text:
#=
module Gaugefields_4D_mpi_module
    using LinearAlgebra
    import ..AbstractGaugefields_module:AbstractGaugefields,Shifted_Gaugefields,shift_U,
                        Adjoint_Gaugefields,set_wing_U!,Abstractfields,construct_staple!,clear_U!,
                        calculate_Plaquette
    import Base
    import ..Gaugefields_4D_module:Gaugefields_4D

    using MPI
    =#

    const comm = MPI.COMM_WORLD

    """
    `Gaugefields_4D_wing_mpi{NC} <: Gaugefields_4D{NC}`

    MPI version of SU(N) Gauge fields in four dimensional lattice.
    """
    struct Gaugefields_4D_wing_mpi{NC} <: Gaugefields_4D{NC}
        U::Array{ComplexF64,6}
        NX::Int64
        NY::Int64
        NZ::Int64
        NT::Int64
        NDW::Int64
        NV::Int64
        NC::Int64
        PEs::NTuple{4,Int64}
        PN::NTuple{4,Int64}
        mpiinit::Bool
        myrank::Int64
        nprocs::Int64
        myrank_xyzt::NTuple{4,Int64}
        mpi::Bool
        verbose_print::Verbose_print

        function Gaugefields_4D_wing_mpi(NC::T,NDW::T,NX::T,NY::T,NZ::T,NT::T,PEs;mpiinit=true,
                                                            verbose_level = 2) where T<: Integer
            NV = NX*NY*NZ*NT
            @assert NX % PEs[1] == 0 "NX % PEs[1] should be 0. Now NX = $NX and PEs = $PEs"
            @assert NY % PEs[2] == 0 "NY % PEs[2] should be 0. Now NY = $NY and PEs = $PEs"
            @assert NZ % PEs[3] == 0 "NZ % PEs[3] should be 0. Now NZ = $NZ and PEs = $PEs"
            @assert NT % PEs[4] == 0 "NT % PEs[4] should be 0. Now NT = $NT and PEs = $PEs"

            PN = (NX ÷ PEs[1],
                    NY ÷ PEs[2],
                    NZ ÷ PEs[3],
                    NT ÷ PEs[4],
            )

            if mpiinit == false
                MPI.Init()
                mpiinit = true
            end

            comm = MPI.COMM_WORLD

            nprocs = MPI.Comm_size(comm)
            @assert prod(PEs) == nprocs "num. of MPI process should be prod(PEs). Now nprocs = $nprocs and PEs = $PEs"
            myrank = MPI.Comm_rank(comm)

            verbose_print = Verbose_print(verbose_level,myid = myrank)

            myrank_xyzt = get_myrank_xyzt(myrank,PEs)

            #println("Hello world, I am $(MPI.Comm_rank(comm)) of $(MPI.Comm_size(comm))")

            U = zeros(ComplexF64,NC,NC,PN[1]+2NDW,PN[2]+2NDW,PN[3]+2NDW,PN[4]+2NDW)
            #U = Array{Array{ComplexF64,6}}(undef,4)
            #for μ=1:4
            #    U[μ] = zeros(ComplexF64,NC,NC,NX+2NDW,NY+2NDW,NZ+2NDW,NT+2NDW)
            #end
            mpi = true
            return new{NC}(U,NX,NY,NZ,NT,NDW,NV,NC,Tuple(PEs),PN,mpiinit,myrank,nprocs,myrank_xyzt,mpi,verbose_print)
        end
    end

    function get_myrank_xyzt(myrank,PEs)
        #myrank = (((myrank_t)*PEs[3]+myrank_z)*PEs[2] + myrank_y)*PEs[1] + myrank_x
        myrank_x = myrank % PEs[1]
        i = (myrank - myrank_x) ÷ PEs[1]
        myrank_y = i % PEs[2]
        i = (i - myrank_y) ÷ PEs[2]
        myrank_z = i % PEs[3]
        myrank_t = (i - myrank_z) ÷ PEs[3]

        return myrank_x,myrank_y,myrank_z,myrank_t
    end

    function get_myrank(myrank_xyzt,PEs)
        @inbounds return (((myrank_xyzt[4])*PEs[3]+myrank_xyzt[3])*PEs[2] + myrank_xyzt[2])*PEs[1] + myrank_xyzt[1]
    end

    function get_myrank(U::T) where T <: Gaugefields_4D_wing_mpi
        return U.myrank
    end

    function get_myrank(U::Array{T,1}) where T <: Gaugefields_4D_wing_mpi
        return U[1].myrank
    end

    function get_nprocs(U::T) where T <: Gaugefields_4D_wing_mpi
        return  U.nprocs
    end

    function get_nprocs(U::Array{T,1}) where T <: Gaugefields_4D_wing_mpi
        return U[1].nprocs
    end

    function calc_rank_and_indices(x::Gaugefields_4D_wing_mpi,ix,iy,iz,it)
        pex = (ix-1) ÷ x.PN[1]
        ix_local = (ix-1) % x.PN[1] + 1

        pey = (iy-1) ÷ x.PN[2]
        iy_local = (iy-1) % x.PN[2] + 1

        pez = (iz-1) ÷ x.PN[3]
        iz_local = (iz-1) % x.PN[3] + 1

        pet = (it-1) ÷ x.PN[4]
        it_local = (it-1) % x.PN[4] + 1
        myrank = get_myrank((pex,pey,pez,pet),x.PEs)
        return myrank,ix_local,iy_local,iz_local,it_local
    end

    function barrier(x::Gaugefields_4D_wing_mpi)
        MPI.Barrier(comm)
    end

    function Base.setindex!(x::Gaugefields_4D_wing_mpi,v,i1,i2,i3,i4,i5,i6)
        error("Each element can not be accessed by global index in $(typeof(x)). Use setvalue! function")
        #x.U[i1,i2,i3 + x.NDW,i4 + x.NDW,i5 + x.NDW,i6 + x.NDW] = v
    end

    function Base.getindex(x::Gaugefields_4D_wing_mpi,i1,i2,i3,i4,i5,i6)
        error("Each element can not be accessed by global index in $(typeof(x)) Use getvalue function")
        #return x.U[i1,i2,i3 .+ x.NDW,i4 .+ x.NDW,i5 .+ x.NDW,i6 .+ x.NDW]
    end

    function Base.setindex!(x::Adjoint_Gaugefields{T},v,i1,i2,i3,i4,i5,i6) where T <: Gaugefields_4D_wing_mpi #U'
        error("type $(typeof(U)) has no setindex method. This type is read only.")
        #x.U[i1,i2,i3 + x.NDW,i4 + x.NDW,i5 + x.NDW,i6 + x.NDW] = v
    end

    function Base.getindex(x::Adjoint_Gaugefields{T},i1,i2,i3,i4,i5,i6) where T <: Gaugefields_4D_wing_mpi #U'
        error("Each element can not be accessed by global index in $(typeof(x)) Use getvalue function")
        #return x.U[i1,i2,i3 .+ x.NDW,i4 .+ x.NDW,i5 .+ x.NDW,i6 .+ x.NDW]
    end

    #=
    function Base.setindex!(U::Shifted_Gaugefields{T,4},v,i1,i2,i3,i4,i5,i6)  where T <: Gaugefields_4D_wing_mpi
        error("type $(typeof(U)) has no setindex method. This type is read only.")
    end

    function Base.getindex(U::Shifted_Gaugefields{T,4},i1,i2,i3,i4,i5,i6) where T <: Gaugefields_4D_wing_mpi
        error("Each element can not be accessed by global index in $(typeof(x)) Use getvalue function")
    end
    =#

    #=

    function Base.getindex(U::Adjoint_Gaugefields{Shifted_Gaugefields{T,4}},i1,i2,i3,i4,i5,i6) where T <: Gaugefields_4D_wing_mpi
        error("Each element can not be accessed by global index in $(typeof(x)) Use getvalue function")
    end

    function Base.setindex!(U::Adjoint_Gaugefields{Shifted_Gaugefields{T,4}},v,i1,i2,i3,i4,i5,i6)  where T <: Gaugefields_4D_wing_mpi
        error("type $(typeof(U)) has no setindex method. This type is read only.")
    end
    =#

    @inline function getvalue(x::Gaugefields_4D_wing_mpi,i1,i2,i3,i4,i5,i6)
        @inbounds return x.U[i1,i2,i3 .+ x.NDW,i4 .+ x.NDW,i5 .+ x.NDW,i6 .+ x.NDW]
    end

    @inline  function setvalue!(x::Gaugefields_4D_wing_mpi,v,i1,i2,i3,i4,i5,i6)
        @inbounds  x.U[i1,i2,i3 + x.NDW,i4 + x.NDW,i5 + x.NDW,i6 + x.NDW] = v
    end

    @inline  function getvalue(x::Adjoint_Gaugefields{T},i1,i2,i3,i4,i5,i6) where T <: Abstractfields
        @inbounds  return conj(getvalue(x.parent,i2,i1,i3,i4,i5,i6))
    end

    @inline  function setvalue!(x::Adjoint_Gaugefields{T},v,i1,i2,i3,i4,i5,i6) where T <: Abstractfields
        error("type $(typeof(x)) has no setindex method. This type is read only.")
    end

    @inline  function getvalue(U::Shifted_Gaugefields{T,4},i1,i2,i3,i4,i5,i6) where T <: Gaugefields_4D_wing_mpi
        @inbounds  return U.parent.U[i1,i2,i3 .+ U.parent.NDW .+ U.shift[1],i4 .+ U.parent.NDW .+ U.shift[2],i5 .+ U.parent.NDW .+ U.shift[3],i6 .+ U.parent.NDW .+ U.shift[4]]
    end

    @inline  function setvalue!(U::Shifted_Gaugefields{T,4},v,i1,i2,i3,i4,i5,i6) where T <: Gaugefields_4D_wing_mpi
        error("type $(typeof(U)) has no setindex method. This type is read only.")
    end

    #=

    @inline  function getvalue(U::Adjoint_Gaugefields{Shifted_Gaugefields{T,4}},i1,i2,i3,i4,i5,i6) where T <: Gaugefields_4D_wing_mpi
        return conj(getvalue(U.parent,i2,i1,i3,i4,i5,i6))
    end

    @inline function setvalue!(U::Adjoint_Gaugefields{Shifted_Gaugefields{T,4}},v,i1,i2,i3,i4,i5,i6) where T <: Gaugefields_4D_wing_mpi
        error("type $(typeof(U)) has no setindex method. This type is read only.")
    end

    =#

    function identityGaugefields_4D_wing_mpi(NC,NX,NY,NZ,NT,NDW,PEs;mpiinit = true,verbose_level = 2,randomnumber="Random")
        U = Gaugefields_4D_wing_mpi(NC,NDW,NX,NY,NZ,NT,PEs,mpiinit = mpiinit,verbose_level = verbose_level)
        v = 1

        for it=1:U.PN[4]
            for iz=1:U.PN[3]
                for iy=1:U.PN[2]
                    for ix=1:U.PN[1]
                        @simd for ic=1:NC
                            setvalue!(U,v,ic,ic,ix,iy,iz,it)
                        end
                    end
                end
            end
        end
        #println("setwing")
        set_wing_U!(U)

        return U
    end

    function randomGaugefields_4D_wing_mpi(NC,NX,NY,NZ,NT,NDW,PEs;mpiinit = true,verbose_level= 2,randomnumber="Random")
        U = Gaugefields_4D_wing_mpi(NC,NDW,NX,NY,NZ,NT,PEs,mpiinit = mpiinit,verbose_level = verbose_level)
        v = 1

        for it=1:U.PN[4]
            for iz=1:U.PN[3]
                for iy=1:U.PN[2]
                    for ix=1:U.PN[1]
                        for jc=1:NC
                            @simd for ic=1:NC
                                v = rand()-0.5 + im*(rand()-0.5)
                                setvalue!(U,v,ic,jc,ix,iy,iz,it)
                            end
                        end
                    end
                end
            end
        end
        #println("setwing")
        normalize_U!(U)
        set_wing_U!(U)

        return U
    end

    function clear_U!(U::Gaugefields_4D_wing_mpi{NC}) where NC
        for it=1:U.PN[4]
            for iz=1:U.PN[3]
                for iy=1:U.PN[2]
                    for ix=1:U.PN[1]
                        for jc=1:NC
                            @simd for ic=1:NC
                                v = 0
                                @inbounds setvalue!(U,v,ic,jc,ix,iy,iz,it)
                                #@inbounds Uμ[k1,k2,ix,iy,iz,it] = 0
                            end
                        end
                    end
                end
            end
        end
        set_wing_U!(U)
    end

    function clear_U!(U::Gaugefields_4D_wing_mpi{NC},iseven::Bool) where NC
        for it=1:U.PN[4]
            for iz=1:U.PN[3]
                for iy=1:U.PN[2]
                    for ix=1:U.PN[1]
                        evenodd = ifelse( (ix+iy+iz+it) % 2 ==0, true,false)
                        if evenodd == iseven
                            for k2=1:NC
                                for k1=1:NC
                                    v = 0
                                    @inbounds setvalue!(U,v,k1,k2,ix,iy,iz,it)
                                end
                            end
                        end
                    end
                end
            end
        end
        set_wing_U!(U)
    end

    function clear_U!(U::Gaugefields_4D_wing_mpi{NC},filternumber::N,filterindex::N) where {NC,N <: Integer}
        for it=1:U.PN[4]
            for iz=1:U.PN[3]
                for iy=1:U.PN[2]
                    for ix=1:U.PN[1]
                        filter = ((ix+iy+iz+it)) % filternumber
                        #evenodd = ifelse( (ix+iy+iz+it) % filternumber ==0, true,false)
                        if filter == filterindex
                            for k2=1:NC
                                for k1=1:NC
                                    v = 0
                                    @inbounds setvalue!(U,v,k1,k2,ix,iy,iz,it)
                                end
                            end
                        end
                    end
                end
            end
        end
        set_wing_U!(U)
    end

    function add_U!(c::Gaugefields_4D_wing_mpi{NC},a::T1) where {NC,T1 <: Abstractfields}
        for it=1:c.PN[4]
            for iz=1:c.PN[3]
                for iy=1:c.PN[2]
                    for ix=1:c.PN[1]

                            for k2=1:NC
                                @simd for k1=1:NC
                                    av = getvalue(a,k1,k2,ix,iy,iz,it)
                                    cv = getvalue(c,k1,k2,ix,iy,iz,it)
                                    v = cv + av
                                    setvalue!(c,v,k1,k2,ix,iy,iz,it)
                                    #c[k1,k2,ix,iy,iz,it] += a[k1,k2,ix,iy,iz,it]
                                end
                            end
                    end
                end
            end
        end
        #set_wing_U!(c)
    end

    function add_U!(c::Gaugefields_4D_wing_mpi{NC},a::T1,iseven::Bool) where {NC,T1 <: Abstractfields}
        @inbounds for it=1:c.PN[4]
            for iz=1:c.PN[3]
                for iy=1:c.PN[2]
                    for ix=1:c.PN[1]
                        evenodd = ifelse( (ix+iy+iz+it) % 2 ==0, true,false)
                        if evenodd == iseven
                            for k2=1:NC
                                @simd for k1=1:NC
                                    av = getvalue(a,k1,k2,ix,iy,iz,it)
                                    cv = getvalue(c,k1,k2,ix,iy,iz,it)
                                    v = cv + av
                                    setvalue!(c,v,k1,k2,ix,iy,iz,it)
                                    #c[k1,k2,ix,iy,iz,it] += a[k1,k2,ix,iy,iz,it]
                                end
                            end
                        end
                    end
                end
            end
        end
        #set_wing_U!(c)
    end

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

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

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

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

    function substitute_U!(U::Gaugefields_4D_wing_mpi{NC},b::T2) where {NC, T2 <: Abstractfields}
        for it=1:U.PN[4]
            for iz=1:U.PN[3]
                for iy=1:U.PN[2]
                    for ix=1:U.PN[1]
                        for k2=1:NC
                            for k1=1:NC
                                v = getvalue(b,k1,k2,ix,iy,iz,it)
                                #v = b[k1,k2,ix,iy,iz,it]
                                @inbounds setvalue!(U,v,k1,k2,ix,iy,iz,it)
                            end
                        end
                    end
                end
            end
        end
        set_wing_U!(U)

    end

    function substitute_U!(U::Gaugefields_4D_wing_mpi{NC},b::T2,iseven::Bool) where {NC, T2 <: Abstractfields}
        for it=1:U.PN[4]
            for iz=1:U.PN[3]
                for iy=1:U.PN[2]
                    for ix=1:U.PN[1]
                        evenodd = ifelse( (ix+iy+iz+it) % 2 ==0, true,false)
                        if evenodd == iseven
                            for k2=1:NC
                                for k1=1:NC
                                    v = getvalue(b,k1,k2,ix,iy,iz,it)
                                    #v = b[k1,k2,ix,iy,iz,it]
                                    @inbounds setvalue!(U,v,k1,k2,ix,iy,iz,it)
                                end
                            end
                        end
                    end
                end
            end
        end
        set_wing_U!(U)

    end

    function map_U!(U::Gaugefields_4D_wing_mpi{NC},f!::Function,V::Gaugefields_4D_wing_mpi{NC},iseven::Bool) where {NC}

        A = zeros(ComplexF64,NC,NC)
        B = zeros(ComplexF64,NC,NC)
        for it=1:U.PN[4]
            for iz=1:U.PN[3]
                for iy=1:U.PN[2]
                    for ix=1:U.PN[1]
                        evenodd = ifelse( (ix+iy+iz+it) % 2 ==0, true,false)
                        if evenodd == iseven
                            for k2=1:NC
                                for k1=1:NC

                                    A[k1,k2] = getvalue(V,k1,k2,ix,iy,iz,it)
                                    B[k1,k2] = getvalue(U,k1,k2,ix,iy,iz,it)
                                end
                            end
                            f!(B,A)
                            for k2=1:NC
                                for k1=1:NC
                                    v = B[k1,k2]
                                    setvalue!(U,v,k1,k2,ix,iy,iz,it)
                                    #U[k1,k2,ix,iy,iz,it] = B[k1,k2]
                                end
                            end
                        end
                    end
                end
            end
        end
        set_wing_U!(U)
    end

    function map_U_sequential!(U::Gaugefields_4D_wing_mpi{NC},f!::Function,Uin) where {NC}
        error("The function map_U_sequential! can not be used with MPI")
    end

    struct Shifted_Gaugefields_4D_mpi{NC,outside} <: Shifted_Gaugefields{NC,4}
        parent::Gaugefields_4D_wing_mpi{NC}
        #parent::T
        shift::NTuple{4,Int8}
        NX::Int64
        NY::Int64
        NZ::Int64
        NT::Int64
        NDW::Int64

        #function Shifted_Gaugefields(U::T,shift,Dim) where {T <: AbstractGaugefields}
        function Shifted_Gaugefields_4D_mpi(U::Gaugefields_4D_wing_mpi{NC},shift) where NC
            outside = check_outside(U.NDW,shift)
            return new{NC,outside}(U,shift,U.NX,U.NY,U.NZ,U.NT,U.NDW)
        end
    end

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

        return Shifted_Gaugefields_4D_mpi(U,shift)
    end

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

    @inline function getvalue(U::Shifted_Gaugefields_4D_mpi{NC,false},i1,i2,i3,i4,i5,i6) where NC
        #function Base.getindex(U::Shifted_Gaugefields{T,4},i1,i2,i3,i4,i5,i6) where T <: Gaugefields_4D_wing
        @inbounds return getvalue(U.parent,i1,i2,i3 .+ U.shift[1],i4 .+ U.shift[2],i5 .+ U.shift[3],i6 .+ U.shift[4])
    end

    @inline function getvalue(U::Shifted_Gaugefields_4D_mpi{NC,true},i1,i2,i3,i4,i5,i6) where NC
        i3_new = i3 + U.shift[1]
        i3_new += ifelse(i3_new > U.NX + U.NDW,-U.NX,0)
        i3_new += ifelse(i3_new < 1 - U.NDW,U.NX,0)
        i4_new = i4 + U.shift[2]
        i4_new += ifelse(i4_new > U.NY + U.NDW,-U.NY,0)
        i4_new += ifelse(i4_new < 1 - U.NDW,U.NY,0)
        i5_new = i5 + U.shift[3]
        i5_new += ifelse(i5_new > U.NZ + U.NDW,-U.NZ,0)
        i5_new += ifelse(i5_new < 1 - U.NDW,U.NZ,0)
        i6_new = i6 + U.shift[4]
        i6_new += ifelse(i6_new > U.NT + U.NDW,-U.NT,0)
        i6_new += ifelse(i6_new < 1 - U.NDW,U.NT,0)
        #function Base.getindex(U::Shifted_Gaugefields{T,4},i1,i2,i3,i4,i5,i6) where T <: Gaugefields_4D_wing
        @inbounds return getvalue(U.parent,i1,i2,i3_new ,i4_new ,i5_new ,i6_new )
    end

    function normalize_U!(U::Gaugefields_4D_wing_mpi{NC}) where NC

        A = zeros(ComplexF64,NC,NC)

        for it=1:U.PN[4]
            for iz=1:U.PN[3]
                for iy=1:U.PN[2]
                    for ix=1:U.PN[1]
                        for jc=1:NC
                            @simd for ic=1:NC
                                A[ic,jc] = getvalue(U,ic,jc,ix,iy,iz,it)
                            end
                        end
                        gramschmidt!(A)

                        for jc=1:NC
                            @simd for ic=1:NC
                                v = A[ic,jc]
                                setvalue!(U,v,ic,jc,ix,iy,iz,it)
                            end
                        end
                    end
                end
            end
        end
        set_wing_U!(U)

    end

    function Base.similar(U::T) where T <: Gaugefields_4D_wing_mpi
        Uout = Gaugefields_4D_wing_mpi(U.NC,U.NDW,U.NX,U.NY,U.NZ,U.NT,U.PEs,mpiinit=U.mpiinit,verbose_level = U.verbose_print.level)
        #identityGaugefields_4D_wing(U.NC,U.NX,U.NY,U.NZ,U.NT,U.NDW)
        return Uout
    end

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

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

        s = 0
        for it=1:PN[4]
            for iz=1:PN[3]
                for iy=1:PN[2]
                    for ix=1:PN[1]
                        @simd for k=1:NC
                            s += getvalue(a,k,k,ix,iy,iz,it)
                            #println(a[k,k,ix,iy,iz,it])
                        end
                    end
                end
            end
        end

        s = MPI.Allreduce(s,MPI.SUM,comm)

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

    function partial_tr(a::Gaugefields_4D_wing_mpi{NC},μ) where NC
        error("Polyakov loop is not supported with MPI yet.")
        PN =a.PN

        if μ == 1
            s = 0
            ix = 1
            for it=1:PN[4]
                for iz=1:PN[3]
                    for iy=1:PN[2]
                        #for ix=1:NX
                            @simd for k=1:NC
                                s += getvalue(a,k,k,ix,iy,iz,it)
                                #println(a[k,k,ix,iy,iz,it])
                            end

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

        end

        s = MPI.Allreduce(s,MPI.SUM,comm)

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

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

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

    function LinearAlgebra.mul!(c::Gaugefields_4D_wing_mpi{NC},a::T1,b::T2,iseven::Bool) where {NC,T1 <: Abstractfields,T2 <: Abstractfields}
        @assert NC != 2 && NC != 3 "This function is for NC != 2,3"
        NT = c.NT
        NZ = c.NZ
        NY = c.NY
        NX = c.NX
        PN = c.PN
        for it=1:PN[4]
            for iz=1:PN[3]
                for iy=1:PN[2]
                    for ix=1:PN[1]
                        evenodd = ifelse( (ix+iy+iz+it) % 2 ==0, true,false)
                        if evenodd == iseven

                            for k2=1:NC
                                for k1=1:NC
                                    v = 0
                                    setvalue!(c,v,k1,k2,ix,iy,iz,it)
                                    #c[k1,k2,ix,iy,iz,it] = 0

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

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

                                @simd for k3=1:NC
                                    av = getvalue(a,k1,k3,ix,iy,iz,it)
                                    bv = getvalue(b,k3,k2,ix,iy,iz,it)
                                    cv = getvalue(c,k1,k2,ix,iy,iz,it)

                                    v = cv + av + bv
                                    setvalue!(c,v,k1,k2,ix,iy,iz,it)
                                    #c[k1,k2,ix,iy,iz,it] += a[k1,k3,ix,iy,iz,it]*b[k3,k2,ix,iy,iz,it]
                                end
                            end
                        end
                    end
                end
            end
        #end
        set_wing_U!(c)
    end

    function LinearAlgebra.mul!(c::Gaugefields_4D_wing_mpi{3},a::T1,b::T2) where {NC,T1 <: Abstractfields,T2 <: Abstractfields}
        #@assert NC != 2 && NC != 3 "This function is for NC != 2,3"
        NT = c.NT
        NZ = c.NZ
        NY = c.NY
        NX = c.NX
        PN = c.PN
        for it=1:PN[4]
            for iz=1:PN[3]
                for iy=1:PN[2]
                    for ix=1:PN[1]
                        a11 = getvalue(a,1,1,ix,iy,iz,it)
                        a21 = getvalue(a,2,1,ix,iy,iz,it)
                        a31 = getvalue(a,3,1,ix,iy,iz,it)
                        a12 = getvalue(a,1,2,ix,iy,iz,it)
                        a22 = getvalue(a,2,2,ix,iy,iz,it)
                        a32 = getvalue(a,3,2,ix,iy,iz,it)
                        a13 = getvalue(a,1,3,ix,iy,iz,it)
                        a23 = getvalue(a,2,3,ix,iy,iz,it)
                        a33 = getvalue(a,3,3,ix,iy,iz,it)
                        b11 = getvalue(b,1,1,ix,iy,iz,it)
                        b21 = getvalue(b,2,1,ix,iy,iz,it)
                        b31 = getvalue(b,3,1,ix,iy,iz,it)
                        b12 = getvalue(b,1,2,ix,iy,iz,it)
                        b22 = getvalue(b,2,2,ix,iy,iz,it)
                        b32 = getvalue(b,3,2,ix,iy,iz,it)
                        b13 = getvalue(b,1,3,ix,iy,iz,it)
                        b23 = getvalue(b,2,3,ix,iy,iz,it)
                        b33 = getvalue(b,3,3,ix,iy,iz,it)

                        v = (a11*b11+a12*b21+a13*b31)
                        setvalue!(c,v,1,1,ix,iy,iz,it)
                        v =  (a21*b11+a22*b21+a23*b31)
                        setvalue!(c,v,2,1,ix,iy,iz,it)
                        v = (a31*b11+a32*b21+a33*b31)
                        setvalue!(c,v,3,1,ix,iy,iz,it)
                        v = (a11*b12+a12*b22+a13*b32)
                        setvalue!(c,v,1,2,ix,iy,iz,it)
                        v = (a21*b12+a22*b22+a23*b32)
                        setvalue!(c,v,2,2,ix,iy,iz,it)
                        v = (a31*b12+a32*b22+a33*b32)
                        setvalue!(c,v,3,2,ix,iy,iz,it)
                        v = (a11*b13+a12*b23+a13*b33)
                        setvalue!(c,v,1,3,ix,iy,iz,it)
                        v =  (a21*b13+a22*b23+a23*b33)
                        setvalue!(c,v,2,3,ix,iy,iz,it)
                        v = (a31*b13+a32*b23+a33*b33)
                        setvalue!(c,v,3,3,ix,iy,iz,it)
                    end
                end
            end
        end
        #set_wing_U!(c)
    end

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

                            a11 = getvalue(a,1,1,ix,iy,iz,it)
                            a21 = getvalue(a,2,1,ix,iy,iz,it)
                            a31 = getvalue(a,3,1,ix,iy,iz,it)
                            a12 = getvalue(a,1,2,ix,iy,iz,it)
                            a22 = getvalue(a,2,2,ix,iy,iz,it)
                            a32 = getvalue(a,3,2,ix,iy,iz,it)
                            a13 = getvalue(a,1,3,ix,iy,iz,it)
                            a23 = getvalue(a,2,3,ix,iy,iz,it)
                            a33 = getvalue(a,3,3,ix,iy,iz,it)
                            b11 = getvalue(b,1,1,ix,iy,iz,it)
                            b21 = getvalue(b,2,1,ix,iy,iz,it)
                            b31 = getvalue(b,3,1,ix,iy,iz,it)
                            b12 = getvalue(b,1,2,ix,iy,iz,it)
                            b22 = getvalue(b,2,2,ix,iy,iz,it)
                            b32 = getvalue(b,3,2,ix,iy,iz,it)
                            b13 = getvalue(b,1,3,ix,iy,iz,it)
                            b23 = getvalue(b,2,3,ix,iy,iz,it)
                            b33 = getvalue(b,3,3,ix,iy,iz,it)

                            v = (a11*b11+a12*b21+a13*b31)
                            setvalue!(c,v,1,1,ix,iy,iz,it)
                            v =  (a21*b11+a22*b21+a23*b31)
                            setvalue!(c,v,2,1,ix,iy,iz,it)
                            v = (a31*b11+a32*b21+a33*b31)
                            setvalue!(c,v,3,1,ix,iy,iz,it)
                            v = (a11*b12+a12*b22+a13*b32)
                            setvalue!(c,v,1,2,ix,iy,iz,it)
                            v = (a21*b12+a22*b22+a23*b32)
                            setvalue!(c,v,2,2,ix,iy,iz,it)
                            v = (a31*b12+a32*b22+a33*b32)
                            setvalue!(c,v,3,2,ix,iy,iz,it)
                            v = (a11*b13+a12*b23+a13*b33)
                            setvalue!(c,v,1,3,ix,iy,iz,it)
                            v =  (a21*b13+a22*b23+a23*b33)
                            setvalue!(c,v,2,3,ix,iy,iz,it)
                            v = (a31*b13+a32*b23+a33*b33)
                            setvalue!(c,v,3,3,ix,iy,iz,it)
                        end
                    end
                end
            end
        end
        #set_wing_U!(c)
    end

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

                        a12 = getvalue(a,1,2,ix,iy,iz,it)
                        a22 = getvalue(a,2,2,ix,iy,iz,it)

                        b11 = getvalue(b,1,1,ix,iy,iz,it)
                        b21 = getvalue(b,2,1,ix,iy,iz,it)

                        b12 = getvalue(b,1,2,ix,iy,iz,it)
                        b22 = getvalue(b,2,2,ix,iy,iz,it)

                        v = a11*b11+a12*b21
                        setvalue!(c,v,1,1,ix,iy,iz,it)
                        v = a21*b11+a22*b21
                        setvalue!(c,v,2,1,ix,iy,iz,it)

                        v = a11*b12+a12*b22
                        setvalue!(c,v,1,2,ix,iy,iz,it)
                        v = a21*b12+a22*b22
                        setvalue!(c,v,2,2,ix,iy,iz,it)
                        #v = a31*b12+a32*b22

                    end
                end
            end
        end
        #set_wing_U!(c)
    end

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

                            a11 = getvalue(a,1,1,ix,iy,iz,it)
                            a21 = getvalue(a,2,1,ix,iy,iz,it)

                            a12 = getvalue(a,1,2,ix,iy,iz,it)
                            a22 = getvalue(a,2,2,ix,iy,iz,it)

                            b11 = getvalue(b,1,1,ix,iy,iz,it)
                            b21 = getvalue(b,2,1,ix,iy,iz,it)

                            b12 = getvalue(b,1,2,ix,iy,iz,it)
                            b22 = getvalue(b,2,2,ix,iy,iz,it)

                            v = a11*b11+a12*b21
                            setvalue!(c,v,1,1,ix,iy,iz,it)
                            v = a21*b11+a22*b21
                            setvalue!(c,v,2,1,ix,iy,iz,it)

                            v = a11*b12+a12*b22
                            setvalue!(c,v,1,2,ix,iy,iz,it)
                            v = a21*b12+a22*b22
                            setvalue!(c,v,2,2,ix,iy,iz,it)
                            v = a31*b12+a32*b22
                        end

                    end
                end
            end
        end
        #set_wing_U!(c)
    end

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

    function LinearAlgebra.mul!(c::Gaugefields_4D_wing_mpi{2},a::T1,b::T2,α::Ta,β::Tb) where {NC,T1 <: Abstractfields,T2 <: Abstractfields,Ta <: Number, Tb <: Number}
        #@assert NC != 2 && NC != 3 "This function is for NC != 2,3"
        NT = c.NT
        NZ = c.NZ
        NY = c.NY
        NX = c.NX
        PN = c.PN
        if β == zero(β)
            if α == one(α)
                mul!(c,a,b)
                return
            end
        end

        @inbounds for it=1:PN[4]
            for iz=1:PN[3]
                for iy=1:PN[2]
                    @simd for ix=1:PN[1]
                        a11 = getvalue(a,1,1,ix,iy,iz,it)
                        a21 = getvalue(a,2,1,ix,iy,iz,it)
                        a12 = getvalue(a,1,2,ix,iy,iz,it)
                        a22 = getvalue(a,2,2,ix,iy,iz,it)

                        b11 = getvalue(b,1,1,ix,iy,iz,it)
                        b21 = getvalue(b,2,1,ix,iy,iz,it)
                        b12 = getvalue(b,1,2,ix,iy,iz,it)
                        b22 = getvalue(b,2,2,ix,iy,iz,it)

                        v = (a11*b11+a12*b21)*α + β*getvalue(c,1,1,ix,iy,iz,it)
                        setvalue!(c,v,1,1,ix,iy,iz,it)
                        v =  (a21*b11+a22*b21)*α + β*getvalue(c,2,1,ix,iy,iz,it)
                        setvalue!(c,v,2,1,ix,iy,iz,it)
                        v = (a11*b12+a12*b22)*α + β*getvalue(c,1,2,ix,iy,iz,it)
                        setvalue!(c,v,1,2,ix,iy,iz,it)
                        v = (a21*b12+a22*b22)*α + β*getvalue(c,2,2,ix,iy,iz,it)
                        setvalue!(c,v,2,2,ix,iy,iz,it)

                    end
                end
            end
        end
        #set_wing_U!(c)
    end

    function LinearAlgebra.mul!(c::Gaugefields_4D_wing_mpi{3},a::T1,b::T2,α::Ta,β::Tb) where {NC,T1 <: Abstractfields,T2 <: Abstractfields,Ta <: Number, Tb <: Number}
        #@assert NC != 2 && NC != 3 "This function is for NC != 2,3"
        NT = c.NT
        NZ = c.NZ
        NY = c.NY
        NX = c.NX
        PN = c.PN
        if β == zero(β)
            if α == one(α)
                mul!(c,a,b)
                return
            end
        end

        @inbounds for it=1:PN[4]
            for iz=1:PN[3]
                for iy=1:PN[2]
                    @simd for ix=1:PN[1]
                        a11 = getvalue(a,1,1,ix,iy,iz,it)
                        a21 = getvalue(a,2,1,ix,iy,iz,it)
                        a31 = getvalue(a,3,1,ix,iy,iz,it)
                        a12 = getvalue(a,1,2,ix,iy,iz,it)
                        a22 = getvalue(a,2,2,ix,iy,iz,it)
                        a32 = getvalue(a,3,2,ix,iy,iz,it)
                        a13 = getvalue(a,1,3,ix,iy,iz,it)
                        a23 = getvalue(a,2,3,ix,iy,iz,it)
                        a33 = getvalue(a,3,3,ix,iy,iz,it)
                        b11 = getvalue(b,1,1,ix,iy,iz,it)
                        b21 = getvalue(b,2,1,ix,iy,iz,it)
                        b31 = getvalue(b,3,1,ix,iy,iz,it)
                        b12 = getvalue(b,1,2,ix,iy,iz,it)
                        b22 = getvalue(b,2,2,ix,iy,iz,it)
                        b32 = getvalue(b,3,2,ix,iy,iz,it)
                        b13 = getvalue(b,1,3,ix,iy,iz,it)
                        b23 = getvalue(b,2,3,ix,iy,iz,it)
                        b33 = getvalue(b,3,3,ix,iy,iz,it)

                        v = (a11*b11+a12*b21+a13*b31)*α + β*getvalue(c,1,1,ix,iy,iz,it)
                        setvalue!(c,v,1,1,ix,iy,iz,it)
                        v =  (a21*b11+a22*b21+a23*b31)*α + β*getvalue(c,2,1,ix,iy,iz,it)
                        setvalue!(c,v,2,1,ix,iy,iz,it)
                        v = (a31*b11+a32*b21+a33*b31)*α + β*getvalue(c,3,1,ix,iy,iz,it)
                        setvalue!(c,v,3,1,ix,iy,iz,it)
                        v = (a11*b12+a12*b22+a13*b32)*α + β*getvalue(c,1,2,ix,iy,iz,it)
                        setvalue!(c,v,1,2,ix,iy,iz,it)
                        v = (a21*b12+a22*b22+a23*b32)*α + β*getvalue(c,2,2,ix,iy,iz,it)
                        setvalue!(c,v,2,2,ix,iy,iz,it)
                        v = (a31*b12+a32*b22+a33*b32)*α + β*getvalue(c,3,2,ix,iy,iz,it)
                        setvalue!(c,v,3,2,ix,iy,iz,it)
                        v = (a11*b13+a12*b23+a13*b33)*α + β*getvalue(c,1,3,ix,iy,iz,it)
                        setvalue!(c,v,1,3,ix,iy,iz,it)
                        v =  (a21*b13+a22*b23+a23*b33)*α + β*getvalue(c,2,3,ix,iy,iz,it)
                        setvalue!(c,v,2,3,ix,iy,iz,it)
                        v = (a31*b13+a32*b23+a33*b33)*α + β*getvalue(c,3,3,ix,iy,iz,it)
                        setvalue!(c,v,3,3,ix,iy,iz,it)

                    end
                end
            end
        end
    end

    function set_wing_U!(u::Array{Gaugefields_4D_wing_mpi{NC},1}) where NC
        for μ=1:4
            set_wing_U!(u[μ])
        end
    end

    function set_wing_U!(u::Gaugefields_4D_wing_mpi{NC}) where NC
        NT = u.NT
        NY = u.NY
        NZ = u.NZ
        NX = u.NX
        NDW = u.NDW
        PEs = u.PEs
        PN = u.PN
        myrank = u.myrank
        myrank_xyzt = u.myrank_xyzt
        myrank_xyzt_send = u.myrank_xyzt

        #X direction
        #Now we send data
        #from NX to 1
        N = PN[2]*PN[3]*PN[4]*NDW*NC*NC
        send_mesg1 = Array{ComplexF64}(undef, N)
        recv_mesg1 = Array{ComplexF64}(undef, N)

        count = 0
        for it=1:PN[4]
            for iz=1:PN[3]
                for iy=1:PN[2]
                    for id=1:NDW
                        for k2=1:NC
                            for k1=1:NC
                                count += 1
                                send_mesg1[count] = getvalue(u,k1,k2,PN[1]+(id-NDW),iy,iz,it)
                                #u[k1,k2,-NDW+id,iy,iz,it] = u[k1,k2,NX+(id-NDW),iy,iz,it]
                            end
                        end
                    end
                end
            end
        end

        px = myrank_xyzt[1] + 1
        px += ifelse(px >= PEs[1],-PEs[1],0)
        myrank_xyzt_send = (px,myrank_xyzt[2],myrank_xyzt[3],myrank_xyzt[4])
        myrank_send1 = get_myrank(myrank_xyzt_send,PEs)
        #=
        for ip=0:u.nprocs-1
            if ip == u.myrank
                println("rank = $myrank, myrank_send1 = $(myrank_send1)")
            end
            MPI.Barrier(comm)

        end
        =#

        sreq1 = MPI.Isend(send_mesg1, myrank_send1, myrank_send1+32, comm) #from left to right 0 -> 1

        N = PN[2]*PN[3]*PN[4]*NDW*NC*NC
        send_mesg2 = Array{ComplexF64}(undef, N)
        recv_mesg2 = Array{ComplexF64}(undef, N)

        count = 0
        for it=1:PN[4]
            for iz=1:PN[3]
                for iy=1:PN[2]
                    for id=1:NDW
                        for k2=1:NC
                            for k1=1:NC
                                count += 1
                                send_mesg2[count] = getvalue(u,k1,k2,id,iy,iz,it)
                            end
                        end
                    end
                end
            end
        end
        px = myrank_xyzt[1] - 1
        px += ifelse(px < 0,PEs[1],0)
        #println("px = $px")
        myrank_xyzt_send = (px,myrank_xyzt[2],myrank_xyzt[3],myrank_xyzt[4])
        myrank_send2 = get_myrank(myrank_xyzt_send,PEs)
        #=
        for ip=0:u.nprocs-1
            if ip == u.myrank
                println("rank = $myrank, myrank_send2 = $(myrank_send2)")
            end
            MPI.Barrier(comm)

        end
        =#

        sreq2 = MPI.Isend(send_mesg2, myrank_send2, myrank_send2+64, comm) #from right to left 0 -> -1

        #=
        myrank = 1: myrank_send1 = 2, myrank_send2 = 0
            sreq1: from 1 to 2 2
            sreq2: from 1 to 0 2
        myrank = 2: myrank_send1 = 3, myrank_send2 = 1
            sreq1: from 2 to 3 3
            sreq2: from 2 to 1 1
            rreq1: from 1 to 2 2 -> sreq1 at myrank 1
            rreq2: from 3 to 2 2
        myrank = 3: myrank_send1 = 4, myrank_send2 = 2
            sreq1: from 3 to 4 4
            sreq2: from 3 to 2 2
        =#

        rreq1 = MPI.Irecv!(recv_mesg1, myrank_send2, myrank+32, comm) #from -1 to 0
        rreq2 = MPI.Irecv!(recv_mesg2, myrank_send1, myrank+64, comm) #from 1 to 0

        stats = MPI.Waitall!([rreq1, sreq1,rreq2,sreq2])
        MPI.Barrier(comm)

        count = 0
        for it=1:PN[4]
            for iz=1:PN[3]
                for iy=1:PN[2]
                    for id=1:NDW
                        for k2=1:NC
                            for k1=1:NC
                                count += 1
                                v = recv_mesg1[count]
                                setvalue!(u,v,k1,k2,-NDW+id,iy,iz,it)
                                #send_mesg1[count] = getvalue(u,k1,k2,PN[1]+(id-NDW),iy,iz,it)
                                #u[k1,k2,-NDW+id,iy,iz,it] = u[k1,k2,NX+(id-NDW),iy,iz,it]
                            end
                        end
                    end
                end
            end
        end

        count = 0
        for it=1:PN[4]
            for iz=1:PN[3]
                for iy=1:PN[2]
                    for id=1:NDW
                        for k2=1:NC
                            for k1=1:NC
                                count += 1
                                v = recv_mesg2[count]
                                setvalue!(u,v,k1,k2,PN[1]+id,iy,iz,it)
                                #u[k1,k2,NX+id,iy,iz,it] = u[k1,k2,id,iy,iz,it]
                                #send_mesg2[count] = getvalue(u,k1,k2,id,iy,iz,it)
                            end
                        end
                    end
                end
            end
        end

        #N = PN[1]*PN[3]*PN[4]*NDW*NC*NC
        N = PN[4]*PN[3]*length(-NDW+1:PN[1]+NDW)*NDW*NC*NC
        send_mesg1 = Array{ComplexF64}(undef, N)
        recv_mesg1 = Array{ComplexF64}(undef, N)
        send_mesg2 = Array{ComplexF64}(undef, N)
        recv_mesg2 = Array{ComplexF64}(undef, N)

        #Y direction
        #Now we send data
        count = 0
        for it=1:PN[4]
            for iz=1:PN[3]
                for ix=-NDW+1:PN[1]+NDW
                    for id=1:NDW
                        for k1=1:NC
                            for k2=1:NC
                                count += 1
                                send_mesg1[count] = getvalue(u,k1,k2,ix,PN[2]+(id-NDW),iz,it)
                                #u[k1,k2,ix,-NDW+id,iz,it] = u[k1,k2,ix,NY+(id-NDW),iz,it]
                            end
                        end
                    end
                end
            end
        end

        py = myrank_xyzt[2] + 1
        py += ifelse(py >= PEs[2],-PEs[2],0)
        myrank_xyzt_send = (myrank_xyzt[1],py,myrank_xyzt[3],myrank_xyzt[4])
        myrank_send1 = get_myrank(myrank_xyzt_send,PEs)
        #println("rank = $rank, myrank_send1 = $(myrank_send1)")
        sreq1 = MPI.Isend(send_mesg1, myrank_send1, myrank_send1+32, comm) #from left to right 0 -> 1

        count = 0
        for it=1:PN[4]
            for iz=1:PN[3]
                for ix=-NDW+1:PN[1]+NDW
                    for id=1:NDW
                        for k1=1:NC
                            for k2=1:NC
                                count += 1
                                send_mesg2[count] = getvalue(u,k1,k2,ix,id,iz,it)
                                #u[k1,k2,ix,NY+id,iz,it] = u[k1,k2,ix,id,iz,it]
                            end
                        end
                    end
                end
            end
        end

        py = myrank_xyzt[2] - 1
        py += ifelse(py < 0,PEs[2],0)
        #println("py = $py")
        myrank_xyzt_send = (myrank_xyzt[1],py,myrank_xyzt[3],myrank_xyzt[4])
        myrank_send2 = get_myrank(myrank_xyzt_send,PEs)
        #println("rank = $rank, myrank_send2 = $(myrank_send2)")
        sreq2 = MPI.Isend(send_mesg2, myrank_send2, myrank_send2+64, comm) #from right to left 0 -> -1

        rreq1 = MPI.Irecv!(recv_mesg1, myrank_send2, myrank+32, comm) #from -1 to 0
        rreq2 = MPI.Irecv!(recv_mesg2, myrank_send1, myrank+64, comm) #from 1 to 0

        stats = MPI.Waitall!([rreq1, sreq1,rreq2,sreq2])

        count = 0
        for it=1:PN[4]
            for iz=1:PN[3]
                for ix=-NDW+1:PN[1]+NDW
                    for id=1:NDW
                        for k1=1:NC
                            for k2=1:NC
                                count += 1
                                v = recv_mesg1[count]
                                setvalue!(u,v,k1,k2,ix,-NDW+id,iz,it)
                                #send_mesg1[count] = getvalue(u,k1,k2,ix,PN[2]+(id-NDW),iz,it)
                                #u[k1,k2,ix,-NDW+id,iz,it] = u[k1,k2,ix,NY+(id-NDW),iz,it]
                            end
                        end
                    end
                end
            end
        end

        count = 0
        for it=1:PN[4]
            for iz=1:PN[3]
                for ix=-NDW+1:PN[1]+NDW
                    for id=1:NDW
                        for k1=1:NC
                            for k2=1:NC
                                count += 1
                                v = recv_mesg2[count]
                                setvalue!(u,v,k1,k2,ix,PN[2]+id,iz,it)
                                #send_mesg2[count] = getvalue(u,k1,k2,ix,id,iz,it)
                                #u[k1,k2,ix,NY+id,iz,it] = u[k1,k2,ix,id,iz,it]
                            end
                        end
                    end
                end
            end
        end

        MPI.Barrier(comm)

        #Z direction
        #Now we send data

        N = NDW*PN[4]*length(-NDW+1:PN[2]+NDW)*length(-NDW+1:PN[1]+NDW)*NC*NC
        send_mesg1 = Array{ComplexF64}(undef, N)
        recv_mesg1 = Array{ComplexF64}(undef, N)
        send_mesg2 = Array{ComplexF64}(undef, N)
        recv_mesg2 = Array{ComplexF64}(undef, N)

        count = 0
        for id=1:NDW
            for it=1:PN[4]
                for iy=-NDW+1:PN[2]+NDW
                    for ix=-NDW+1:PN[1]+NDW
                        for k1=1:NC
                            for k2=1:NC
                                count += 1
                                send_mesg1[count] = getvalue(u,k1,k2,ix,iy,PN[3]+(id-NDW),it)
                                send_mesg2[count] = getvalue(u,k1,k2,ix,iy,id,it)
                                #u[k1,k2,ix,iy,id-NDW,it] = u[k1,k2,ix,iy,NZ+(id-NDW),it]
                                #u[k1,k2,ix,iy,NZ+id,it] = u[k1,k2,ix,iy,id,it]
                            end
                        end
                    end
                end
            end
        end

        pz = myrank_xyzt[3] + 1
        pz += ifelse(pz >= PEs[3],-PEs[3],0)
        myrank_xyzt_send = (myrank_xyzt[1],myrank_xyzt[2],pz,myrank_xyzt[4])
        myrank_send1 = get_myrank(myrank_xyzt_send,PEs)
        #println("rank = $rank, myrank_send1 = $(myrank_send1)")
        sreq1 = MPI.Isend(send_mesg1, myrank_send1, myrank_send1+32, comm) #from left to right 0 -> 1

        pz = myrank_xyzt[3] - 1
        pz += ifelse(pz < 0,PEs[3],0)
        #println("pz = $pz")
        myrank_xyzt_send = (myrank_xyzt[1],myrank_xyzt[2],pz,myrank_xyzt[4])
        myrank_send2 = get_myrank(myrank_xyzt_send,PEs)
        #println("rank = $rank, myrank_send2 = $(myrank_send2)")
        sreq2 = MPI.Isend(send_mesg2, myrank_send2, myrank_send2+64, comm) #from right to left 0 -> -1

        rreq1 = MPI.Irecv!(recv_mesg1, myrank_send2, myrank+32, comm) #from -1 to 0
        rreq2 = MPI.Irecv!(recv_mesg2, myrank_send1, myrank+64, comm) #from 1 to 0

        stats = MPI.Waitall!([rreq1, sreq1,rreq2,sreq2])

        count = 0
        for id=1:NDW
            for it=1:PN[4]
                for iy=-NDW+1:PN[2]+NDW
                    for ix=-NDW+1:PN[1]+NDW
                        for k1=1:NC
                            for k2=1:NC
                                count += 1
                                v = recv_mesg1[count]
                                setvalue!(u,v,k1,k2,ix,iy,id-NDW,it)
                                v = recv_mesg2[count]
                                setvalue!(u,v,k1,k2,ix,iy,PN[3]+id,it)
                                #u[k1,k2,ix,iy,id-NDW,it] = u[k1,k2,ix,iy,NZ+(id-NDW),it]
                                #u[k1,k2,ix,iy,NZ+id,it] = u[k1,k2,ix,iy,id,it]
                            end
                        end
                    end
                end
            end
        end

        MPI.Barrier(comm)

        #T direction
        #Now we send data

        N = NDW*length(-NDW+1:PN[3]+NDW)*length(-NDW+1:PN[2]+NDW)*length(-NDW+1:PN[1]+NDW)*NC*NC
        send_mesg1 = Array{ComplexF64}(undef, N)
        recv_mesg1 = Array{ComplexF64}(undef, N)
        send_mesg2 = Array{ComplexF64}(undef, N)
        recv_mesg2 = Array{ComplexF64}(undef, N)

        count = 0
        for id=1:NDW
            for iz=-NDW+1:PN[3]+NDW
                for iy=-NDW+1:PN[2]+NDW
                    for ix=-NDW+1:PN[1]+NDW
                        for k1=1:NC
                            for k2=1:NC
                                count += 1
                                send_mesg1[count] = getvalue(u,k1,k2,ix,iy,iz,PN[4]+(id-NDW))
                                send_mesg2[count] = getvalue(u,k1,k2,ix,iy,iz,id)
                                #u[k1,k2,ix,iy,iz,id-NDW] = u[k1,k2,ix,iy,iz,PN[4]+(id-NDW)]
                                #u[k1,k2,ix,iy,iz,PN[4]+id] = u[k1,k2,ix,iy,iz,id]
                            end
                        end
                    end
                end
            end
        end

        pt = myrank_xyzt[4] + 1
        pt += ifelse(pt >= PEs[4],-PEs[4],0)
        myrank_xyzt_send = (myrank_xyzt[1],myrank_xyzt[2],myrank_xyzt[3],pt)
        myrank_send1 = get_myrank(myrank_xyzt_send,PEs)
        #println("rank = $rank, myrank_send1 = $(myrank_send1)")
        sreq1 = MPI.Isend(send_mesg1, myrank_send1, myrank_send1+32, comm) #from left to right 0 -> 1

        pt = myrank_xyzt[4] - 1
        pt += ifelse(pt < 0,PEs[4],0)
        #println("pt = $pt")
        myrank_xyzt_send = (myrank_xyzt[1],myrank_xyzt[2],myrank_xyzt[3],pt)
        myrank_send2 = get_myrank(myrank_xyzt_send,PEs)
        #println("rank = $rank, myrank_send2 = $(myrank_send2)")
        sreq2 = MPI.Isend(send_mesg2, myrank_send2, myrank_send2+64, comm) #from right to left 0 -> -1

        rreq1 = MPI.Irecv!(recv_mesg1, myrank_send2, myrank+32, comm) #from -1 to 0
        rreq2 = MPI.Irecv!(recv_mesg2, myrank_send1, myrank+64, comm) #from 1 to 0

        stats = MPI.Waitall!([rreq1, sreq1,rreq2,sreq2])

        count = 0
        for id=1:NDW
            for iz=-NDW+1:PN[3]+NDW
                for iy=-NDW+1:PN[2]+NDW
                    for ix=-NDW+1:PN[1]+NDW
                        for k1=1:NC
                            for k2=1:NC
                                count += 1
                                v = recv_mesg1[count]
                                setvalue!(u,v,k1,k2,ix,iy,iz,id-NDW)
                                v = recv_mesg2[count]
                                setvalue!(u,v,k1,k2,ix,iy,iz,PN[4]+id)

                                #send_mesg1[count] = getvalue(u,k1,k2,ix,iy,iz,PN[4]+(id-NDW))
                                #send_mesg2[count] = getvalue(u,k1,k2,ix,iy,iz,id)
                                #u[k1,k2,ix,iy,iz,id-NDW] = u[k1,k2,ix,iy,iz,PN[4]+(id-NDW)]
                                #u[k1,k2,ix,iy,iz,PN[4]+id] = u[k1,k2,ix,iy,iz,id]
                            end
                        end
                    end
                end
            end
        end
        #error("rr22r")

        MPI.Barrier(comm)

        return
    end

#end


text:
module TestStructuredConcurrency
include("../../../examples/structured_concurrency.jl")
end  # module

ext: jl


text:
@testset "fhir-to-julia.jl" begin
    @testset "fhir_to_julia" begin
        @test FHIRClient.fhir_to_julia(:end) == :end_fhir
        @test FHIRClient.fhir_to_julia(:function) == :function_fhir
        @test FHIRClient.fhir_to_julia(:global) == :global_fhir
        @test FHIRClient.fhir_to_julia(:import) == :import_fhir
        @test FHIRClient.fhir_to_julia(:for) == :for_fhir
        @test FHIRClient.fhir_to_julia(:birthDate) == :birthDate
    end

    @testset "julia_to_fhir" begin
        @test FHIRClient.julia_to_fhir(:end_fhir) == :end
        @test FHIRClient.julia_to_fhir(:function_fhir) == :function
        @test FHIRClient.julia_to_fhir(:global_fhir) == :global
        @test FHIRClient.julia_to_fhir(:import_fhir) == :import
        @test FHIRClient.julia_to_fhir(:for_fhir) == :for
        @test FHIRClient.julia_to_fhir(:birthDate) == :birthDate
    end
end


text:
using RobotZoo
using RobotDynamics
using StaticArrays
using BenchmarkTools
using FiniteDiff
using LinearAlgebra
using Test
using Random
Random.seed!(1)
const RD = RobotDynamics

# Continuous Linear model
n,m = 10,5
A,B = RobotZoo.RandomLinearModels.gencontrollable(n,m, :continuous)
d = randn(10)
model = RobotZoo.LinearModel(A,B,d)
@test RD.dims(model) == (n,m,n)
x,u = rand(model)
@test length(x) == n
@test length(u) == m
@test RD.dynamics(model, x, u) ≈ A*x + B*u + d
xn = similar(x)
RD.dynamics!(model, xn, x, u)
@test xn ≈ A*x + B*u + d

@test RobotZoo.iscontrollable(model)

J = zeros(n, n+m)
z = RD.KnotPoint(x, u, 0, NaN)
RD.jacobian!(RD.InPlace(), RD.UserDefined(), model, J, xn, z)
@test J[:,1:n] ≈ A
@test J[:,n+1:end] ≈ B

Jfd = similar(J)
FiniteDiff.finite_difference_jacobian!(
    Jfd,
    (y,z)->RD.dynamics!(model, y, z[1:n], z[n+1:end]),
    Vector(z.z)
)
@test Jfd ≈ J rtol=1e-6

# Linearized Continuous Model
Random.seed!(1)
model0 = RobotZoo.Cartpole()
xe = [0,pi,0,0]
ue = [0.]
model = RobotZoo.LinearModel(RD.InPlace(), RD.ForwardAD(), model0, xe, ue)
n,m = RD.dims(model)
@test n == 4
@test m == 1
dx = randn(4)
du = randn(1)
@test RD.dynamics(model, dx*0, du*0) ≈ RD.dynamics(model0, xe, ue)
@test norm(RD.dynamics(model, dx*1e-3, du*1e-3) - RD.dynamics(model0, xe, ue)) < 0.1

Jfd = zeros(4,5)
y = zeros(4)
z = RD.KnotPoint{n,m}(xe, ue, 0, NaN)
FiniteDiff.finite_difference_jacobian!(
    Jfd,
    (y,z)->RD.dynamics!(model0, y, z[1:n], z[n+1:end]),
    Vector(z.z)
)
@test Jfd[:,1:n] ≈ RobotZoo.getstatecoeffs(model)
@test Jfd[:,n+1:end] ≈ RobotZoo.getcontrolcoeffs(model)
@test RobotZoo.isstable(model) == false

# Linearize about stable equilibrium
xe = [0,0.0,0,0]
ue = [0.]
model_stable = RobotZoo.LinearModel(RD.InPlace(), RD.ForwardAD(), model0, xe, ue)
@test RobotZoo.isstable(model_stable)

# Use random generator on linear model
n,m = 12,6
model_rand = rand(RobotZoo.LinearModel, n, m)
@test RD.dims(model_rand) == (n,m,n)
@test RobotZoo.iscontrollable(model_rand)

## Discrete Linear Model
n,m = 10,5
h = 0.1
A,B = RobotZoo.RandomLinearModels.gencontrollable(n, m, :discrete)
d = randn(n)
model = RobotZoo.DiscreteLinearModel(h,A,B)
@test RobotZoo.getstatecoeffs(model) == A
@test RobotZoo.getcontrolcoeffs(model) == B
@test RobotZoo.getaffinecoeffs(model) == zeros(n)
@test RobotZoo.iscontrollable(model)

model = RobotZoo.DiscreteLinearModel(h,A,B,d)
@test RobotZoo.getstatecoeffs(model) == A
@test RobotZoo.getcontrolcoeffs(model) == B
@test RobotZoo.getaffinecoeffs(model) == d
@test RobotZoo.iscontrollable(model)

C = randn(n,n)
model = RobotZoo.DiscreteLinearModel(h,A,B,C,d)
@test RobotZoo.getstatecoeffs(model) == A
@test RobotZoo.getcontrolcoeffs(model) == B
@test RobotZoo.getnextstatecoeffs(model) == C
@test RobotZoo.getaffinecoeffs(model) == d

@test_nowarn RobotZoo.checktimestep(model, h)
@test_nowarn RobotZoo.checktimestep(model, h + eps(h))
@test_throws ArgumentError RobotZoo.checktimestep(model, h*1.1)

@test RD.dims(model) == (n,m,n)
x,u = rand(model)
z = RD.KnotPoint(x,u,0.,h)

# Explicit
model_explicit = RobotZoo.DiscreteLinearModel(h,A,B,d)
@test RobotZoo.getnextstatecoeffs(model_explicit) ≈ -I(n)
@test RD.discrete_dynamics(model_explicit, x, u, 0, h) ≈ A*x + B*u + d

J = zeros(n,n+m)
xn = zeros(n)
RD.jacobian!(RD.InPlace(), RD.UserDefined(), model_explicit, J, xn, z)
@test J[:,1:n] ≈ A
@test J[:,n+1:end] ≈ B

# Implicit
model_implicit = RobotZoo.DiscreteLinearModel(h,A,B,C,d)
@test RobotZoo.getnextstatecoeffs(model_implicit) ≈ C
@test RD.discrete_dynamics(model_implicit, x, u, 0, h) ≈ -C\(A*x + B*u + d)

RD.jacobian!(RD.InPlace(), RD.UserDefined(), model_implicit, J, xn, z)
@test J[:,1:n] ≈ -C\A
@test J[:,n+1:end] ≈ -C\B

## Discretize and Linearize
model0 = RobotZoo.Cartpole()
n,m = RD.dims(model0)

# Implicit
dmodel_explicit = RD.DiscretizedDynamics{RD.RK3}(model0)
xe = [0,pi,0,0]
ue = [0.]
model_explicit = RobotZoo.DiscreteLinearModel(
    RD.InPlace(),
    RD.ForwardAD(),
    dmodel_explicit,
    xe, ue, 0., h
)

J = zeros(n, n+m)
z = RD.KnotPoint{n,m}([xe;ue], 0.0, h)
xn = zeros(n)
RD.jacobian!(RD.StaticReturn(), RD.ForwardAD(), dmodel_explicit, J, xn, z)
@test J[:,1:n] ≈ RobotZoo.getstatecoeffs(model_explicit)
@test J[:,n+1:end] ≈ RobotZoo.getcontrolcoeffs(model_explicit)
@test RD.discrete_dynamics(dmodel_explicit, z) ≈
    RobotZoo.getaffinecoeffs(model_explicit)
@test RobotZoo.getnextstatecoeffs(model_explicit) == -I(n) * one(eltype(A))
@test RobotZoo.isstable(model_explicit) == false
@test RobotZoo.iscontrollable(model_explicit)

xe_stable = zeros(n)
model_explicit_stable = RobotZoo.DiscreteLinearModel(
    RD.InPlace(),
    RD.ForwardAD(),
    dmodel_explicit,
    xe_stable, ue, 0., h
)
@test RobotZoo.isstable(model_explicit_stable) == true

# Implicit
dmodel_implicit = RD.DiscretizedDynamics{RD.ImplicitMidpoint}(model0)
model_implicit= RobotZoo.DiscreteLinearModel(
    RD.InPlace(),
    RD.ForwardAD(),
    dmodel_implicit,
    xe, ue, 0., h
)
@test RD.dims(model_implicit) == (n,m,n)

Jc = zeros(n,n+m)
RD.jacobian!(RD.StaticReturn(), RD.ForwardAD(), model0, Jc, xn, z)
A = Jc[:,1:n]
B = Jc[:,n+1:end]
@test RobotZoo.getstatecoeffs(model_implicit) ≈ h/2*A + I
@test RobotZoo.getnextstatecoeffs(model_implicit) ≈ h/2*A - I
@test RobotZoo.getcontrolcoeffs(model_implicit) ≈ h*B

# Exponential integrator
model_linear = RobotZoo.LinearModel(RD.InPlace(), RD.ForwardAD(),
    model0, xe, ue, affine=false)
model_exp = RobotZoo.DiscreteLinearModel(model_linear, h)

# Compare the linear models
dx = randn(n) * 1e-2
du = randn(m) * 1e-2
xn_exp = RD.discrete_dynamics(model_exp, dx, du, 0.0, h)
xn_implicit = RD.discrete_dynamics(model_implicit, dx, du, 0, h) - model_implicit.d
xn_explicit = RD.discrete_dynamics(model_explicit, dx, du, 0, h) - model_explicit.d
@test norm(xn_exp - xn_implicit) < 1e-3
@test norm(xn_exp - xn_explicit) < 1e-3


text:
@testset "module evaluation" begin

function mock_output(str, mod = nothing)
    result_doc = mock_run(str; mod = mod)
    return result_doc.chunks[1].output
end

str = """
```julia
@__MODULE__
```
"""

# in sandbox
@test occursin(r"\#+WeaveSandBox[\#\d]+", mock_output(str))

# in Main
@test strip(mock_output(str, Main)) == "Main"

end # @testset "module evaluation"

@testset "clear_module!" begin

ary = rand(1000000)
size = Base.summarysize(ary)

# simple case
m = Core.eval(@__MODULE__, :(module $(gensym(:WeaveTestModule)) end))
Core.eval(m, :(a = $ary))
Weave.clear_module!(m)
@test Base.summarysize(m) < size

# recursive case
m = Core.eval(@__MODULE__, :(module $(gensym(:WeaveTestModule)) end))
Core.eval(m, :(
    module $(gensym(:WeaveTestSubModule))
    a = $ary
    end
))
Weave.clear_module!(m)
@test Base.summarysize(m) < size

# doesn't work with constants
m = Core.eval(@__MODULE__, :(module $(gensym(:WeaveTestModule)) end))
Core.eval(m, :(const a = $ary))
Weave.clear_module!(m)
@test_broken Base.summarysize(m) < size

end # @testset "clear_module!"


text:
# TODO: consider linked list to avoid concatenation costs

module PilotSCAs

export SCA
export numStates, numActions
export reward, nextStates

export State, Action

using GridInterpolations, DiscreteMDPs

import DiscreteMDPs.DiscreteMDP
import DiscreteMDPs.reward
import DiscreteMDPs.nextStates
import DiscreteMDPs.numStates
import DiscreteMDPs.numActions

using PilotSCAConst

type SCA <: DiscreteMDP

    nStates::Int64
    nActions::Int64
    actions::Vector{Symbol}
    grid::RectangleGrid

    function SCA()

        grid = RectangleGrid(Xs, Ys, Bearings, Speeds, Speeds, Responses, Responses)

        return new(NStates, NActions, Actions, grid)

    end # function SCA

end # type SCA

type State

    x::Float64
    y::Float64
    bearing::Float64
    speedOwnship::Float64
    speedIntruder::Float64

    clearOfConflict::Bool

    respondingOwnship::Bool
    respondingIntruder::Bool

end # type State

type Action

    ownship::Symbol
    intruder::Symbol

end # type Action

function numStates(mdp::SCA)

    return mdp.nStates

end # function numStates

function numActions(mdp::SCA)

    return mdp.nActions

end # function numActions

function reward(mdp::SCA, istate::Int64, iaction::Int64)

    state = State(
        0.0,  # x
        0.0,  # y
        0.0,  # bearing
        10.0,  # speedOwnship
        10.0,  # speedIntruder
        true,  # clearOfConflict
        false,  # respondingOwnship
        false)  # respondingIntruder

    if istate < mdp.nStates
        state = gridState2state(ind2x(mdp.grid, istate))
    end # if

    action = ind2a(mdp.actions, iaction)

    return reward(mdp, state, action)

end # function reward

function reward(mdp::SCA, state::State, action::Action)

    reward = 0.0

    if action.ownship != :clearOfConflict
        reward -= PenConflict
    end # if

    if action.intruder != :clearOfConflict
        reward -= PenConflict
    end # if

    turnOwnship = getTurnAngle(action.ownship, state.respondingOwnship)
    turnIntruder = getTurnAngle(action.intruder, state.respondingIntruder)
    reward -= PenAction * (turnOwnship^2 + turnIntruder^2)

    if !state.clearOfConflict
        minSepSq = Inf
        for ti = 1:DT / DTI
            minSepSq = min(minSepSq, getSepSq(state))
            state = getNextState(state, action, 0.0, 0.0, DTI)
        end # for ti

        if minSepSq < MinSepSq
            reward -= PenMinSep
        end # if

        reward -= PenCloseness * exp(-minSepSq * InvVar)
    end # if

    return reward

end # function reward

function ind2a(actions::Vector{Symbol}, iaction::Int64)

    iOwnship = iaction % length(actions)
    if iOwnship == 0
        iOwnship = 6
    end # if

    iIntruder = (iaction - iOwnship) / length(actions) + 1

    return Action(actions[iOwnship], actions[iIntruder])

end # function ind2a

# Returns turn angle corresponding to action in degrees.
function getTurnAngle(action::Symbol, responding::Bool)

    if !responding
        return 0.0
    elseif action == :clearOfConflict
        return 0.0
    elseif action == :straight
        return 0.0
    elseif action == :left10
        return 10.0
    elseif action == :right10
        return -10.0
    elseif action == :left20
        return 20.0
    elseif action == :right20
        return -20.0
    else
        throw(ArgumentError("illegal action symbol"))
    end # if

end # function getTurnAngle

function getSepSq(state::State)

    if state.clearOfConflict
        return Inf
    else
        return state.x^2 + state.y^2
    end # if

end # function getSepSq

function getNextState(
        state::State,
        action::Action,
        sigmaTurnOwnship::Float64 = 0.0,
        sigmaTurnIntruder::Float64 = 0.0,
        dt::Float64 = DT)

    newState = State(
        state.x,
        state.y,
        state.bearing,
        state.speedOwnship,
        state.speedIntruder,
        state.clearOfConflict,
        state.respondingOwnship,
        state.respondingIntruder)

    if !state.clearOfConflict

        turnOwnship = deg2rad(getTurnAngle(action.ownship, state.respondingOwnship)) + sigmaTurnOwnship
        turnIntruder = deg2rad(getTurnAngle(action.intruder, state.respondingIntruder)) + sigmaTurnIntruder

        if turnOwnship == 0.0 || turnIntruder == 0.0  # straight line path(s)

            if turnIntruder != 0.0  # ownship straight path

                gtan = G * tan(turnIntruder)
                bearingChange = dt * gtan / state.speedIntruder
                radiusIntruder = abs(state.speedIntruder^2 / gtan)

                newX = state.x + radiusIntruder * sign(bearingChange) * (sin(state.bearing) - sin(state.bearing - bearingChange)) - state.speedOwnship * dt
                newY = state.y + radiusIntruder * sign(bearingChange) * (-cos(state.bearing) + cos(state.bearing - bearingChange))
                newBearing = norm_angle(state.bearing + bearingChange)

                if newX < Xmin || newX > Xmax || newY < Ymin || newY > Ymax
                    newState.clearOfConflict = true
                else
                    newState.x = newX
                    newState.y = newY
                    newState.bearing = newBearing
                end # if

            elseif turnOwnship != 0.0  # intruder straight path

                gtan = G * tan(turnOwnship)
                bearingChange = dt * gtan / state.speedOwnship
                radiusOwnship = abs(state.speedOwnship^2 / gtan)

                x = state.x + state.speedIntruder * dt * cos(state.bearing) - radiusOwnship * sign(bearingChange) * sin(bearingChange)
                y = state.y + state.speedIntruder * dt * sin(state.bearing) - radiusOwnship * sign(bearingChange) * (cos(bearingChange) - 1)

                newX = x * cos(bearingChange) + y * sin(bearingChange)
                newY = -x * sin(bearingChange) + y * cos(bearingChange)
                newBearing = norm_angle(state.bearing - bearingChange)

                if newX < Xmin || newX > Xmax || newY < Ymin || newY > Ymax
                    newState.clearOfConflict = true
                else
                    newState.x = newX
                    newState.y = newY
                    newState.bearing = newBearing
                end # if

            else  # both straight paths

                newX = state.x + state.speedIntruder * dt * cos(state.bearing) - state.speedOwnship * dt
                newY = state.y + state.speedIntruder * dt * sin(state.bearing)
                newBearing = norm_angle(state.bearing)

                if newX < Xmin || newX > Xmax || newY < Ymin || newY > Ymax
                    newState.clearOfConflict = true
                else
                    newState.x = newX
                    newState.y = newY
                    newState.bearing = newBearing
                end # if

            end # if

        else  # both curved paths

            gtanOwnship = G * tan(turnOwnship)
            gtanIntruder = G * tan(turnIntruder)

            bearingChangeOwnship = dt * gtanOwnship / state.speedOwnship
            bearingChangeIntruder = dt * gtanIntruder / state.speedIntruder

            radiusOwnship = abs(state.speedOwnship^2 / gtanOwnship)
            radiusIntruder = abs(state.speedIntruder^2 / gtanIntruder)

            x = state.x + radiusIntruder * sign(bearingChangeIntruder) * (sin(state.bearing) - sin(state.bearing - bearingChangeIntruder))
              - radiusOwnship * sign(bearingChangeOwnship) * sin(bearingChangeOwnship)
            y = state.y + radiusIntruder * sign(bearingChangeIntruder) * (-cos(state.bearing) + cos(state.bearing - bearingChangeIntruder))
              - radiusOwnship * sign(bearingChangeOwnship) * (-1 + cos(bearingChangeOwnship))

            newX = x * cos(bearingChangeOwnship) + y * sin(bearingChangeOwnship)
            newY = -x * sin(bearingChangeOwnship) + y * cos(bearingChangeOwnship)
            newBearing = norm_angle(state.bearing + bearingChangeIntruder - bearingChangeOwnship)

            if newX < Xmin || newX > Xmax || newY < Ymin || newY > Ymax
                newState.clearOfConflict = true
            else
                newState.x = newX
                newState.y = newY
                newState.bearing = newBearing
            end # if

        end # if

    end # if

    return newState

end # function getNextState

function norm_angle(angle::Float64)
    return ((angle % (2 * pi)) + 2 * pi) % (2 * pi)
end # function norm_angle

# Returns next states and associated transition probabilities.
function nextStates(mdp::SCA, istate::Int64, iaction::Int64)

    if istate == mdp.nStates
        return [mdp.nStates], [1.0]
    end # if

    state = gridState2state(ind2x(mdp.grid, istate))
    action = ind2a(mdp.actions, iaction)

    if action.ownship == :clearOfConflict && action.intruder == :clearOfConflict

        return sigmaSample(mdp, state, action)

    elseif action.ownship != :clearOfConflict && action.intruder == :clearOfConflict

        if state.respondingOwnship

            return sigmaSample(mdp, state, action)

        else  # !state.respondingOwnship

            # case 1: ownship doesn't respond
            nonresponsiveIndices, nonresponsiveProbs = sigmaSample(mdp, state, action)

            # case 2: ownship responds
            state.respondingOwnship = true
            responsiveIndices, responsiveProbs = sigmaSample(mdp, state, action)

            return [nonresponsiveIndices, responsiveIndices],
                   [nonresponsiveProbs * NonRespondingProb, responsiveProbs * RespondingProb]

        end # if

    elseif action.ownship == :clearOfConflict && action.intruder != :clearOfConflict

        if state.respondingIntruder

            return sigmaSample(mdp, state, action)

        else  # !state.respondingIntruder

            # case 1: intruder doesn't respond
            nonresponsiveIndices, nonresponsiveProbs = sigmaSample(mdp, state, action)

            # case 2: intruder responds
            state.respondingIntruder = true
            responsiveIndices, responsiveProbs = sigmaSample(mdp, state, action)

            return [nonresponsiveIndices, responsiveIndices],
                   [nonresponsiveProbs * NonRespondingProb, responsiveProbs * RespondingProb]

        end # if

    else  # action.ownship != :clearOfConflict && action.intruder != :clearOfConflict

        if state.respondingOwnship && state.respondingIntruder

            return sigmaSample(mdp, state, action)

        elseif !state.respondingOwnship && state.respondingIntruder

            # case 1: ownship doesn't respond
            nonresponsiveIndices, nonresponsiveProbs = sigmaSample(mdp, state, action)

            # case 2: ownship responds
            state.respondingOwnship = true
            responsiveIndices, responsiveProbs = sigmaSample(mdp, state, action)

            return [nonresponsiveIndices, responsiveIndices],
                   [nonresponsiveProbs * NonRespondingProb, responsiveProbs * RespondingProb]

        elseif state.respondingOwnship && !state.respondingIntruder

            # case 1: intruder doesn't respond
            nonresponsiveIndices, nonresponsiveProbs = sigmaSample(mdp, state, action)

            # case 2: intruder responds
            state.respondingIntruder = true
            responsiveIndices, responsiveProbs = sigmaSample(mdp, state, action)

            return [nonresponsiveIndices, responsiveIndices],
                   [nonresponsiveProbs * NonRespondingProb, responsiveProbs * RespondingProb]

        else  # !state.respondingOwnship && !state.respondingOwnship

            # case 1: both ownship and intruder don't respond
            noresponseIndices, noresponseProbs = sigmaSample(mdp, state, action)

            # case 2: ownship responds but not intruder
            state.respondingOwnship = true
            ownshipResponseIndices, ownshipResponseProbs = sigmaSample(mdp, state, action)

            # case 3: both ownship and intruder respond
            state.respondingIntruder = true
            bothResponseIndices, bothResponseProbs = sigmaSample(mdp, state, action)

            # case 4: intruder responds but not ownship
            state.respondingOwnship = false
            intruderResponseIndices, intruderResponseProbs = sigmaSample(mdp, state, action)

            return [
                noresponseIndices,
                ownshipResponseIndices,
                bothResponseIndices,
                intruderResponseIndices],
                [
                noresponseProbs * NonRespondingProb * NonRespondingProb,
                ownshipResponseProbs * RespondingProb * NonRespondingProb,
                bothResponseProbs * RespondingProb * RespondingProb,
                intruderResponseProbs * NonRespondingProb * RespondingProb]

        end # if

    end # if

end # function nextStates

function sigmaSample(mdp::SCA, state::State, action::Action)

    nominalIndices, nominalProbs = nextStatesSigma(mdp, state, action)
    speedIndices, speedProbs = sigmaSpeed(mdp, state, action)
    bankIndices, bankProbs = sigmaBank(mdp, state, action)

    return [
            nominalIndices,
            speedIndices,
            bankIndices],
        [
            nominalProbs * SigmaWeightNominal,
            speedProbs * SigmaWeightOffNominal,
            bankProbs * SigmaWeightOffNominal]

end # function sigmaSample

function sigmaSpeed(mdp::SCA, state::State, action::Action)

    # negative sigma
    state.speedOwnship -= SigmaSpeed
    negIndicesOwnship, negProbsOwnship = nextStatesSigma(mdp, state, action)

    # positive sigma
    state.speedOwnship += 2 * SigmaSpeed
    posIndicesOwnship, posProbsOwnship = nextStatesSigma(mdp, state, action)

    # restore original
    state.speedOwnship -= SigmaSpeed

    # negative sigma
    state.speedIntruder -= SigmaSpeed
    negIndicesIntruder, negProbsIntruder = nextStatesSigma(mdp, state, action)

    # positive sigma
    state.speedIntruder += 2 * SigmaSpeed
    posIndicesIntruder, posProbsIntruder = nextStatesSigma(mdp, state, action)

    # restore original
    state.speedIntruder -= SigmaSpeed

    return [negIndicesOwnship, posIndicesOwnship, negIndicesIntruder, posIndicesIntruder],
           [negProbsOwnship, posProbsOwnship, negProbsIntruder, posProbsIntruder]

end # function sigmaSpeed

function sigmaBank(mdp::SCA, state::State, action::Action)

    sigmaBankVal = SigmaBank

    if !state.respondingOwnship || action.ownship == :clearOfConflict
        sigmaBankVal = SigmaBankCOC
    end # if

    # negative sigma
    negIndicesOwnship, negProbsOwnship = nextStatesSigma(mdp, state, action, -sigmaBankVal, 0.0)

    # positive sigma
    posIndicesOwnship, posProbsOwnship = nextStatesSigma(mdp, state, action, sigmaBankVal, 0.0)

    if !state.respondingIntruder || action.intruder == :clearOfConflict
        sigmaBankVal = SigmaBankCOC
    end # if

    # negative sigma
    negIndicesIntruder, negProbsIntruder = nextStatesSigma(mdp, state, action, 0.0, -sigmaBankVal)

    # positive sigma
    posIndicesIntruder, posProbsIntruder = nextStatesSigma(mdp, state, action, 0.0, sigmaBankVal)

    return [negIndicesOwnship, posIndicesOwnship, negIndicesIntruder, posIndicesIntruder],
           [negProbsOwnship, posProbsOwnship, negProbsIntruder, posProbsIntruder]

end # function sigmaBank

function nextStatesSigma(
        mdp::SCA,
        state::State,
        action::Action,
        sigmaTurnOwnship::Float64 = 0.0,
        sigmaTurnIntruder::Float64 = 0.0)

    trueNextState = getNextState(state, action, sigmaTurnOwnship, sigmaTurnIntruder)
    gridNextState = state2gridState(trueNextState)

    if trueNextState.clearOfConflict
        return [mdp.nStates], [1.0]
    else
        return interpolants(mdp.grid, gridNextState)
    end # if

end # function nextStatesSigma

function state2gridState(state::State)

    # indicator values for responsiveness: false == 0, true == 1
    respondingOwnship = 0.0
    respondingIntruder = 0.0

    if state.respondingOwnship
        respondingOwnship = 1.0
    end # if

    if state.respondingIntruder
        respondingIntruder = 1.0
    end # if

    return [
        state.x,
        state.y,
        state.bearing,
        state.speedOwnship,
        state.speedIntruder,
        respondingOwnship,
        respondingIntruder]

end # function state2gridState

function index2state(mdp::SCA, stateIndices::Vector{Int64})

    states = Array(State, length(stateIndices))

    for index = 1:length(stateIndices)
        states[index] = gridState2state(ind2x(mdp.grid, index))
    end # for index

    return states

end # function index2state

function gridState2state(gridState::Vector{Float64})

    # indicator values for responsiveness: false == 0, true == 1
    respondingOwnship = false
    respondingIntruder = false

    if gridState[6] == 1.0
        respondingOwnship = true
    end # if

    if gridState[7] == 1.0
        respondingIntruder = true
    end # if

    return State(
        gridState[1],  # x
        gridState[2],  # y
        gridState[3],  # bearing
        gridState[4],  # speedOwnship
        gridState[5],  # speedIntruder
        false,  # clearOfConflict
        respondingOwnship,  # respondingOwnship
        respondingIntruder)  # respondingIntruder

end # function gridState2state

end # module PilotSCAs


text:
"""
rse(
    holo::AbstractMatrix{T},
    windows_size::Tuple{Integer,Integer},
    interval::Tuple{Integer,Integer},
    P::FFTW.FFTWPlan,
    scale::Integer;
    shift = false
    ) -> AbstractMatrix{T} where T <: AbstractFloat

`冗余散斑降噪法`的算法实现，具体可见`doi: 10.1364/AO.390500`

# Arguments
- holo: 全息图（非RGB类型而是AbstractFloat类型）
- windows_size: 掩膜窗口大小
- interval: 窗口偏移量
- holo: 原始全息图
- P: 其类型为FFTW.FFTWPlan，表示重建采用`fft`傅立叶变换
- scale: 由于零级像强度过大，故需要用`scale`拉升±1级像的强度
- shift: 是否旋转图像，默认为否（除非是想展示完整的重建像，否则不建议将`shift`设为`true`，因为其内存和时间开销比较大）
"""
function rse(holo::AbstractMatrix{T}, windows_size::Tuple{Integer,Integer},
    interval::Tuple{Integer,Integer}, P::FFTW.FFTWPlan, scale::Integer; shift = false) where T <: AbstractFloat

    # 计算子全息图的个数并进行预分配内存
    nx, ny = (size(holo) .- windows_size) .÷ interval
    S = (nx + 1) * (ny + 1)
    holo_tensor = zeros(T, size(holo, 1), size(holo, 2), S)

    # 获取子全息图
    make_sub_holo!(holo_tensor, holo, windows_size, interval, (nx, ny))

    if typeof(P) <: FFTW.rFFTWPlan
        re_tensor = similar(holo, size(holo, 1) ÷ 2 + 1, size(holo, 2), S)
        rse_img = zeros(T, size(holo, 1) ÷ 2 + 1, size(holo, 2))
    elseif typeof(P) <: FFTW.cFFTWPlan
        re_tensor = similar(holo, size(holo, 1), size(holo, 2), S)
        rse_img = zeros(T, size(holo, 1), size(holo, 2))
    else
        @error "P的类型不正确..."
    end

    # 对子全息图进行重建
    reconst_tensor!(re_tensor, holo_tensor, P, scale; shift = shift)

    # 求其均值
    rse_core!(rse_img, re_tensor)

    return rse_img

end

"""
    rse_core!(rse_img::AbstractMatrix{T}, re_tensor::AbstractArray{T,3}) -> AbstractMatrix where T <: AbstractFloat

`冗余散斑降噪法`的核心算法，对子再现像序列进行去冗余而后叠加平均，并存进`rse_img`中

# Arguments
- rse_img: 存放结果的矩阵
- re_tensor: 子再现像序列
"""
function rse_core!(rse_img::AbstractMatrix{T}, re_tensor::AbstractArray{T,3}) where T <: AbstractFloat
    dims1, dims2, dims3 = size(re_tensor)
    @inbounds Threads.@threads for ind in CartesianIndices((1:dims1, 1:dims2))
        i, j = ind.I
        hash_tab = Set{T}()
        clock = 0
        @inbounds for k in 1: dims3
            if re_tensor[i, j, k] ∉ hash_tab
                push!(hash_tab, re_tensor[i, j, k])
                rse_img[i, j] += re_tensor[i, j, k]
            clock += 1
            end
        end

        rse_img[i, j] /= clock
    end

    return rse_img
end


text:
N_A = 6.02214076E23
m_e = 5.48579909070E-04
# Comment

max_stars_repo_name: wamcghee/NautilusMC

max_issues_repo_name: wamcghee/NautilusMC

max_forks_repo_name: wamcghee/NautilusMC


text:
function detect(net, meta, img; thresh=0.5, hier_thresh=0.5, nms=0)
    num = Cint(0)
    pnum = Ref(num)
    network_predict_image(net, img)
    dets = get_network_boxes(net, img.w, img.h, thresh, hier_thresh, C_NULL, 0, pnum, 0)
    num = pnum[]
    if nms > 0
        do_nms_sort(dets, num, meta.classes, nms)
    end
    res = []
    detections = unsafe_wrap(Array, dets, num)
    names = unsafe_wrap(Array, meta.names, meta.classes)
    for j in 1:num
        for i in 1:meta.classes
            prob = unsafe_wrap(Array, detections[j].prob, meta.classes)
            if prob[i] > 0
                b = detections[j].bbox
                nameTag = names[i]
                push!(res,(nameTag, prob[i], (b.x, b.y, b.w, b.h)))
            end
        end
    end
    free_detections(dets, num)
    return res
end


text:
using Documenter, ScatteringTransform

makedocs(;
    modules=[ScatteringTransform],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/dsweber2/ScatteringTransform.jl/blob/{commit}{path}#L{line}",
    sitename="ScatteringTransform.jl",
    authors="dsweber2",
    assets=String[],
)

deploydocs(;
    repo="github.com/dsweber2/ScatteringTransform.jl",
)


text:
#! /bin/bash
julia --color=yes make.jl
open build/index.html


text:
using SequentialMonteCarlo
using RNGPool
using SMCExamples.LinearGaussian: LGTheta, Float64Particle, kalmanlogZ,
  defaultLGModel, makeLGModel
using MonteCarloMarkovKernels
using StaticArrays
using Test
using Random
using LinearAlgebra
import Statistics.mean

setRNGs(0)
lgModel, theta, ys, ko = defaultLGModel(10)

const truex0 = theta.x0
const truev0 = theta.v0
const trueC = theta.C

@inline function toLGTheta(v::SVector{3, Float64})
  return LGTheta(v[1], v[2], trueC, v[3], truex0, truev0)
end

t0 = SVector{3, Float64}(theta.A, theta.Q, theta.R)
const sigmaProp = SMatrix{3, 3, Float64}(Matrix{Float64}(I, 3, 3))

@inline function lglogprior(theta::LGTheta)
  if theta.A < 0 || theta.A > 10 return -Inf end
  if theta.Q < 0 || theta.Q > 10 return -Inf end
  if theta.R < 0 || theta.R > 10 return -Inf end
  return 0.0
end

function makelgsmcltd(ys::Vector{Float64}, N::Int64, nthreads::Int64)
  smcio = SMCIO{Float64Particle, Nothing}(N, length(ys), nthreads, false, 0.5)
  function ltd(in::SVector{3, Float64})
    theta::LGTheta = toLGTheta(in)
    lp::Float64 = lglogprior(theta)
    if lp == -Inf return -Inf end
    model::SMCModel = makeLGModel(theta, ys)
    smc!(model, smcio)
    return smcio.logZhats[length(ys)]
  end
end

function makelgkalmanltd(ys::Vector{Float64})
  function ltd(in::SVector{3, Float64})
    theta::LGTheta = toLGTheta(in)
    lp::Float64 = lglogprior(theta)
    if lp == -Inf return -Inf end
    return kalmanlogZ(theta, ys)
  end
end

logtargetSMC = makelgsmcltd(ys, 128, Threads.nthreads())
logtargetKalman = makelgkalmanltd(ys)

PSMC = makeAMKernel(logtargetSMC, sigmaProp)
PKalman = makeAMKernel(logtargetKalman, sigmaProp)

Random.seed!(12345)

chainSMC = simulateChainProgress(PSMC, t0, 1024*32)

chainKalman = simulateChain(PKalman, t0, 1024*1024)

@test mean(chainSMC) ≈ mean(chainKalman) rtol=0.1
@test MonteCarloMarkovKernels.cov(chainSMC) ≈
  MonteCarloMarkovKernels.cov(chainKalman) rtol=0.1


text:
function __init__mutagenesis()
    ORIGINAL_LINK = "https://relational.fit.cvut.cz/dataset/Mutagenesis"
    DATA_LINK = "https://raw.githubusercontent.com/CTUAvastLab/datasets/main/mutagenesis"
    DEPNAME = "Mutagenesis"
    DATA = "data.json"
    METADATA = "meta.json"

    register(DataDep(
        DEPNAME,
        """
        Dataset: The $DEPNAME dataset.
        Website: $ORIGINAL_LINK
        License: CC0
        """,
        "$DATA_LINK/" .* [DATA, METADATA],
        "80ec1716217135e1f2e0b5a61876c65184e2014e64551103c41e174775ca207c"
    ))
end

"""
    Mutagenesis(; split, dir=nothing)

The `Mutagenesis` dataset comprises 188 molecules trialed for mutagenicity on Salmonella typhimurium, available from
 [relational.fit.cvut.cz](https://relational.fit.cvut.cz/dataset/Mutagenesis) and
 [CTUAvastLab/datasets](https://github.com/CTUAvastLab/datasets/tree/main/mutagenesis).

Set `split` to `:train`, `:val`, `:test`, or `:all`, to select the training,
validation, test partition respectively or the whole dataset.
The `indexes` field in the result contains the indexes of the partition in the
full dataset.

Website: https://relational.fit.cvut.cz/dataset/Mutagenesis
License: CC0

```julia-repl
julia> using MLDatasets: Mutagenesis

julia> dataset = Mutagenesis(split=:train)
Mutagenesis dataset:
  split : train
  indexes : 100-element Vector{Int64}
  features : 100-element Vector{Dict{Symbol, Any}}
  targets : 100-element Vector{Int64}

julia> dataset[1].features
Dict{Symbol, Any} with 5 entries:
  :lumo  => -1.246
  :inda  => 0
  :logp  => 4.23
  :ind1  => 1
  :atoms => Dict{Symbol, Any}[Dict(:element=>"c", :bonds=>Dict{Symbol, Any}[Dict(:element=>"c", :bond_type=>7, :charge=>-0.117, :atom_type=>22), Dict(:element=>"h", :bond_type=>1, :charge=>0.142, :atom_type=>3)…

julia> dataset[1].targets
1

julia> dataset = Mutagenesis(split=:all)
Mutagenesis dataset:
  split : all
  indexes : 188-element Vector{Int64}
  features : 188-element Vector{Dict{Symbol, Any}}
  targets : 188-element Vector{Int64}
```
"""
struct Mutagenesis <: SupervisedDataset
    metadata::Dict{Symbol, Any}
    split::Symbol
    indexes::Vector{Int}    # indexes of the split in the full dataset
    features::Vector{Dict{Symbol, Any}}
    targets::Vector{Int}
end

function Mutagenesis(; split::Symbol, dir=nothing)
    DEPNAME = "Mutagenesis"
    DATA = "data.json"
    METADATA = "meta.json"

    @assert split ∈ [:train, :val, :test, :all]

    data_path = datafile(DEPNAME, DATA, dir)
    metadata_path = datafile(DEPNAME, METADATA, dir)
    samples = open(JSON3.read, data_path)
    metadata = open(JSON3.read, metadata_path)
    labelkey = metadata["label"]
    targets = map(i -> i[labelkey], samples)
    features = map(x->delete!(copy(x), Symbol(labelkey)), samples)
    val_num = metadata["val_samples"]
    test_num = metadata["test_samples"]
    train_idxs = 1:length(samples)-val_num-test_num
    val_idxs = length(samples)-val_num-test_num+1:length(samples)-test_num
    test_idxs = length(samples)-test_num+1:length(samples)
    indexes = split == :train ? train_idxs :
              split == :val ? val_idxs :
              split == :test ? test_idxs : 1:length(samples)
    metadata = Dict(metadata)
    metadata[:data_path] = data_path
    metadata[:metadata_path] = metadata_path

    Mutagenesis(metadata, split, indexes, features[indexes], targets[indexes])
end

# deprecated in v0.6
function Base.getproperty(::Type{Mutagenesis}, s::Symbol)
    if s == :traindata
        @warn "Mutagenesis.traindata() is deprecated, use `Mutagenesis(split=:train)` instead."
        return (; dir=nothing) -> begin
                    d = Mutagenesis(; split=:train, dir)
                    d.features, d.targets
                end
    elseif s == :valdata
        @warn "Mutagenesis.valdata() is deprecated, use `Mutagenesis(split=:val)` instead."
        return (; dir=nothing) -> begin
                    d = Mutagenesis(; split=:val, dir)
                    d.features, d.targets
                end
    elseif s == :testdata
        @warn "Mutagenesis.testdata() is deprecated, use `Mutagenesis(split=:test)` instead."
        return (; dir=nothing) -> begin
                    d = Mutagenesis(; split=:test, dir)
                    d.features, d.targets
                end
    else
        return getfield(Mutagenesis, s)
    end
end


text:
using Temporal, Test

@testset "OHLC" begin
    data = TS(rand(252, 8))
    data.fields = [:Open, :High, :Low, :Close, :Volume, :Settle, :AdjClose, :Last]
    @test has_close(data, allow_settle=true,  allow_last=true)
    @test has_close(data, allow_settle=false, allow_last=true)
    @test has_close(data, allow_settle=false, allow_last=false)
    @test has_close(data, allow_settle=true, allow_last=false)
    @test has_open(data)
    @test has_low(data)
    @test has_high(data)
    @test has_volume(data)
    @test op(data).fields[1] == :Open
    @test hi(data).fields[1] == :High
    @test lo(data).fields[1] == :Low
    @test vo(data).fields[1] == :Volume
    @test cl(data).fields[1] == :AdjClose
    @test cl(data[:,1:end-2]).fields[1] == :Settle
    @test cl(data, use_adj=false, allow_settle=false).fields[1] == :Close
    @test ohlc(data) == data[[:Open,:High,:Low,:AdjClose]]
    @test ohlcv(data) == data[[:Open,:High,:Low,:AdjClose,:Volume]]
    @test hlc(data) == data[[:High,:Low,:AdjClose]]
    @test hl(data) == data[[:High,:Low]]
    @test hl2(data).values == apply(hl(data), fun=mean).values
    @test hlc3(data).values == apply(hlc(data), fun=mean).values
    @test ohlc4(data).values == apply(ohlc(data), fun=mean).values
end


text:
__precompile__(true)

module FourRoomsContNNExperiment

using Reproduce
using JLD2
using Resampling
import Resampling.ExpUtils: algorithm_settings!, build_algorithm_dict
import Resampling.ExpUtils.FluxUtil
using Resampling.ExpUtils.FourRoomsContUtil
# import JuliaRL
using Random
using Statistics
import Flux: Descent
import ProgressMeter
import StatsBase
using JuliaRL.FeatureCreators

function mave(preds, truth, args...)
    return mean(abs.(preds .- truth))
end

function get_target_policy_matrix(env::FourRooms, μ::Resampling.AbstractPolicy)
    tmp = zeros(4, 11, 11)
    for state in Resampling.get_states(env)
        tmp[:, state] .= [get(μ, state, action, nothing, nothing, nothing) for action in 1:4]
    end
    return tmp
end

function exp_settings(as::ArgParseSettings = ArgParseSettings(exc_handler=Reproduce.ArgParse.debug_handler, exit_after_help=false))
    @add_arg_table as begin
        # Basic Experiment Arguments
        "--exp_loc"
        arg_type=String
        required=true
        "--seed"
        help = "Starting seed (initialize seed + run)"
        default = 0
        arg_type = Int64
        "--run"
        help = "Run"
        default = 1
        arg_type = Int64
        "--numinter"
        help = "Number of interactions"
        default = 1000
        arg_type = Int64
        "--train_gap"
        help = "Number of interactions in-between training iterations"
        default = 1
        arg_type = Int64
        "--warm_up"
        help = "number of interactions to wait to begin learning"
        default = 0
        arg_type = Int64
        "--buffersize"
        help = "Size of buffer"
        default = 15000
        arg_type = Int64
        "--batchsize"
        help = "Size of batch"
        default = 16
        arg_type = Int64
        "--eval_points"
        help="Number of points to evaluate over"
        default = 50
        arg_type = Int64
        "--eval_steps"
        help="Number of points to evaluate over"
        default = 1
        arg_type = Int64
        "--working"
        action=:store_true
        "--compress"
        action=:store_true
    end

    FluxUtil.opt_settings!(as)
    algorithm_settings!(as)
    env_settings!(as)

    return as
end

function main_experiment(args::Vector{String})

    as = exp_settings()
    parsed = parse_args(args, as)
    if parsed == nothing
        return
    end
    save_loc=""
    if !parsed["working"]
        create_info!(parsed, parsed["exp_loc"]; filter_keys=["verbose", "working", "exp_loc"])
        save_loc = Reproduce.get_save_dir(parsed)
        if isfile(joinpath(save_loc, "results.jld2"))
            return
        end
    end

    # num_runs = parsed["numruns"]
    train_gap = parsed["train_gap"]
    warm_up = parsed["warm_up"]
    num_interactions = parsed["numinter"]
    buffer_size = parsed["buffersize"]
    batch_size = parsed["batchsize"]

    μ = get_policy(parsed)
    gvf = FourRoomsContUtil.GVFS[parsed["gvf"]]
    opt = FluxUtil.get_optimizer(parsed)

    rng = MersenneTwister(parsed["seed"] + parsed["run"])

    env = FourRoomsCont(0.1, 0.001)
    # eval_states = [Resampling.random_start_state(env, rng) for i in 1:parsed["eval_points"]]
    eval_states = FourRoomsContUtil.sample_according_to_dmu(env, parsed["policy"], parsed["eval_points"]; rng=rng)
    eval_rets = [mean(Resampling.MonteCarloReturn(env, gvf, start_state, 100; rng=rng)) for start_state in eval_states]

    agent = NNFourRoomsContAgent(μ, gvf, opt, train_gap, buffer_size, batch_size, warm_up, parsed, size(env);
                                 rng=rng, normalize=(x::Array{Float64, 1})->(x./11))

    # agent = TCFourRoomsContAgent(μ, gvf, 64, 8, α_arr, train_gap, buffer_size, batch_size, warm_up, parsed, size(env))

    error_dict = Dict{String, Array{Float64}}()
    for key in keys(agent.algo_dict)
        error_dict[key] = zeros(Float32, Int64(floor(num_interactions/parsed["eval_steps"])))
    end

    _, s_t = start!(env; rng=rng)
    action = start!(agent, s_t; rng=rng)
    predict_dict = Dict{String, Array{Float64, 1}}()

    eval_step = 1

    # ProgressMeter.@showprogress 0.1 "Step: " for step in 1:num_interactions
    for step in 1:num_interactions

        # Get experience from environment.
        _, s_tp1, r, terminal = step!(env, action; rng=rng)
        action = step!(agent, s_tp1, r, terminal; rng=rng)

        # calculate error
        if (step-1)%parsed["eval_steps"] == 0
            predict!(agent, eval_states, predict_dict)
            for key in keys(agent.algo_dict)
                error_dict[key][eval_step] = Float32(mave(predict_dict[key], eval_rets))
            end
            eval_step += 1
        end
    end
    # # println(value_function)

    if parsed["working"]
        return error_dict
    else
        results = error_dict
        if parsed["compress"]
            JLD2.save(
                JLD2.FileIO.File(JLD2.FileIO.DataFormat{:JLD2},
                                 joinpath(save_loc, "results.jld2")),
                Dict("results"=>results); compress=true)
        else
            @save joinpath(save_loc, "results.jld2") error_dict
        end
    end

end

Base.@ccallable function julia_main(ARGS::Vector{String})::Cint
    main_experiment(ARGS)
    return 0
end

end