text:
import Base

import Gaugefields:comm,setvalue!
import Gaugefields:barrier
import Gaugefields.AbstractGaugefields_module:getvalue

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

    function Data_sent_fermion(N,NC;NG = 4)
        data = zeros(ComplexF64,NC,NG,N)
        count = 0
        positions = zeros(Int64,N)

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

"""
Struct for WilsonFermion
"""
struct WilsonFermion_4D_nowing_mpi{NC} <: WilsonFermion_4D{NC} #AbstractFermionfields_4D{NC}
    f::Array{ComplexF64,6}
    NC::Int64
    NX::Int64
    NY::Int64
    NZ::Int64
    NT::Int64
    NG::Int64
    NDW::Int64
    Dirac_operator::String
    PEs::NTuple{4,Int64} #number of processes in each dimension
    PN::NTuple{4,Int64} #number of sites in each process
    mpiinit::Bool
    myrank::Int64
    nprocs::Int64
    myrank_xyzt::NTuple{4,Int64}
    mpi::Bool
    fshifted::Array{ComplexF64,6}
    tempmatrix::Array{ComplexF64,3}
    positions::Vector{Int64}
    send_ranks::Dict{Int64,Data_sent_fermion{NC}}
    win::MPI.Win
    win_i::MPI.Win
    win_1i::MPI.Win
    countvec::Vector{Int64}
    otherranks::Vector{Int64}
    win_other::MPI.Win
    your_ranks::Matrix{Int64}
    comm::MPI.Comm
    #BoundaryCondition::Vector{Int8}

    function WilsonFermion_4D_nowing_mpi(NC::T,NX::T,NY::T,NZ::T,NT::T,PEs;comm = MPI.COMM_WORLD) where T<: Integer
        NG = 4
        NDW = 0
        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],
            )

        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)

        myrank_xyzt = get_myrank_xyzt(myrank,PEs)

        #@assert NDW == 1 "only NDW = 1 is supported. Now NDW = $NDW"
        f = zeros(ComplexF64,NC,NG,PN[1]+2NDW,PN[2]+2NDW,PN[3]+2NDW,PN[4]+2NDW) #note: ic,ialpha,ix,iy,it,iz
        fshifted = zero(f)

        Dirac_operator = "Wilson"
        mpi = true
        mpiinit = true

        tempmatrix = zeros(ComplexF64,NC,NG,prod(PN))
        positions = zeros(Int64,prod(PN))
        send_ranks = Dict{Int64,Data_sent_fermion{NC}}()
        mpi = true
        win = MPI.Win_create(tempmatrix,comm)
        win_i = MPI.Win_create(positions,comm)
        countvec = zeros(Int64,1)
        win_1i = MPI.Win_create(countvec,comm)

        otherranks = zeros(Int64,nprocs)
        otherranks .= 0
        win_other = MPI.Win_create(otherranks,comm)
        your_ranks = zeros(Int64,nprocs,nprocs)

        return new{NC}(f,NC,NX,NY,NZ,NT,NG,NDW,Dirac_operator,Tuple(PEs),PN,mpiinit,myrank,nprocs,myrank_xyzt,mpi,
                    fshifted,tempmatrix,positions,send_ranks,
                    win,win_i,win_1i,countvec,otherranks,win_other,your_ranks,comm)
    end

end

function get_myrank(x::T) where T <: WilsonFermion_4D_nowing_mpi
    return x.myrank
end

function get_nprocs(x::T) where T <: WilsonFermion_4D_nowing_mpi
    return  x.nprocs
end

function barrier(x::T) where T <: WilsonFermion_4D_nowing_mpi
    MPI.Barrier(x.comm)
end

function Base.similar(x::T) where T <:  WilsonFermion_4D_nowing_mpi
    return WilsonFermion_4D_nowing_mpi(x.NC,x.NX,x.NY,x.NZ,x.NT,x.PEs,comm = x.comm)
end

function Base.setindex!(x::WilsonFermion_4D_nowing_mpi{NC},v,i1,i2,i3,i4,i5,i6)  where {NC}
    error("Each element can not be accessed by global index in $(typeof(x)). Use setvalue! function")

end

function Base.getindex(x::WilsonFermion_4D_nowing_mpi{NC},i1,i2,i3,i4,i5,i6) where {NC}
    @warn "Each element can not be accessed by global index in $(typeof(x)) Use getvalue function"
    return getvalue(x,i1,i2,i3,i4,i5,i6)
    #error("Each element can not be accessed by global index in $(typeof(x)) Use getvalue function")

end

function setindex_global!(x::WilsonFermion_4D_nowing_mpi{NC},v,ic,ix,iy,iz,it,ialpha)  where {NC}
    i1 = ic
    i2 = ialpha
    i3 = ix
    i4 = iy
    i5 = iz
    i6 = it

    PN = x.PN
    PEs = x.PEs
    ii3 = i3 + ifelse(i3 < 1,x.NX,0) + ifelse(i3 > x.NX,-x.NX,0)
    ii4 = i4 + ifelse(i4 < 1,x.NY,0) + ifelse(i4 > x.NY,-x.NY,0)
    ii5 = i5 + ifelse(i5 < 1,x.NZ,0) + ifelse(i5 > x.NZ,-x.NZ,0)
    ii6 = i6 + ifelse(i6 < 1,x.NT,0) + ifelse(i6 > x.NT,-x.NT,0)
    #i = myrank_xyz*PN + i_local
    myrank_x = ii3 ÷ PN[1]
    myrank_y = ii4 ÷ PN[2]
    myrank_z = ii5 ÷ PN[3]
    myrank_t = ii6 ÷ PN[4]
    myrank = (((myrank_t)*PEs[3]+myrank_z)*PEs[2] + myrank_y)*PEs[1] + myrank_x
    if myrank == x.myrank
        ilocal_3 = ((ii3-1) % PN[1]) + 1
        ilocal_4 = ((ii4-1) % PN[2]) + 1
        ilocal_5 = ((ii5-1) % PN[3]) + 1
        ilocal_6 = ((ii6-1) % PN[4]) + 1
        setvalue!(x,v,i1,i2,ilocal_3,ilocal_4,ilocal_5,ilocal_6)
    end
end

@inline function getvalue(x::WilsonFermion_4D_nowing_mpi{NC},i1,i2,i3,i4,i5,i6) where {NC}
    @inbounds return x.f[i1,i2,i3 ,i4 ,i5 ,i6]
end

@inline function  getvalue(F::Adjoint_fermionfields{T},i1,i2,i3,i4,i5,i6) where T <: Abstractfermion  #F'
    @inbounds return conj(getvalue(F.parent,i1,i2,i3,i4,i5,i6))
end

@inline  function setvalue!(x::WilsonFermion_4D_nowing_mpi{NC},v,i1,i2,i3,i4,i5,i6) where {NC}
    @inbounds x.f[i1,i2 ,i3 ,i4 ,i5 ,i6] = v
end

struct Shifted_fermionfields_4D_nowing_mpi{NC,T} <: Shifted_fermionfields{NC,4}
    parent::T
    #parent::T
    shift::NTuple{4,Int8}
    NC::Int64

    #function Shifted_Gaugefields(U::T,shift,Dim) where {T <: AbstractGaugefields}
    function Shifted_fermionfields_4D_nowing_mpi(F::WilsonFermion_4D_nowing_mpi{NC},shift;boundarycondition = boundarycondition_default ) where {NC}
        shifted_fermion!(F,boundarycondition,shift)
        return new{NC,typeof(F)}(F,shift,NC)
    end
end

function getvalue(F::Shifted_fermionfields_4D_nowing_mpi,i1,i2,i3,i4,i5,i6)
    @inbounds return F.parent.fshifted[i1,i2,i3,i4,i5,i6]
end

function getvalue(F::Shifted_fermionfields_4D_nowing_mpi,i1::N,i2::N,i3::N,i4::N,i5::N,i6::N)  where {N <: Integer}
    @inbounds return F.parent.fshifted[i1,i2,i3,i4,i5,i6]
end

function shift_fermion(F::WilsonFermion_4D_nowing_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_fermionfields_4D_nowing_mpi(F,shift)
end

function shift_fermion(F::TF,shift::NTuple{Dim,T}) where {Dim,T <: Integer,TF <: WilsonFermion_4D_nowing_mpi}
    return Shifted_fermionfields_4D_nowing_mpi(F,shift)
end

function update_sent_data!(send_ranks,N,ix,iy,iz,it,ix_shifted,iy_shifted,iz_shifted,it_shifted,
    PEs,myrank_xyzt,xP,yP,zP,tP,x::WilsonFermion_4D_nowing_mpi{NC},factor) where {NC}
    NG = x.NG
    tempmatrix_mini = view(x.tempmatrix,1:NC,1:NG,1)

    px = myrank_xyzt[1] + xP
    while px >= PEs[1]
        px += -PEs[1]
    end
    while px < 0
        px += PEs[1]
    end

    py = myrank_xyzt[2] + yP
    while py >= PEs[2]
        py += -PEs[2]
    end
    while py < 0
        py += PEs[2]
    end

    pz = myrank_xyzt[3] + zP
    while pz >= PEs[3]
        pz += -PEs[3]
    end
    while pz < 0
        pz += PEs[3]
    end

    pt = myrank_xyzt[4] + tP
    while pt >= PEs[4]
        pt += -PEs[4]
    end
    while pt < 0
        pt += PEs[4]
    end

    myrank_xyzt_send = (px,py,pz,pt)

    myrank_send = get_myrank(myrank_xyzt_send,PEs)
    #println("send ",myrank_send)

    for jc=1:NG
        @simd for ic=1:NC
            #v = getvalue(U,ic,jc,ix_shifted_back,iy_shifted_back,iz_shifted_back,it_shifted_back)
            #v = getvalue(U,ic,jc,ix_shifted,iy_shifted,iz_shifted,it_shifted)
            v = factor*getvalue(x,ic,jc,ix,iy,iz,it)
            tempmatrix_mini[ic,jc] = v
        end
    end
    #disp = ((((it-1)*x.PN[3] + iz-1)*x.PN[2] + iy-1)*x.PN[1] + ix-1)*NC*NC
    #disp = ((((it_shifted-1)*x.PN[3] + iz_shifted-1)*x.PN[2] + iy_shifted-1)*x.PN[1] + ix_shifted-1)*NC*NC
    #println(myrank_send)
    disp = (((it_shifted-1)*x.PN[3] + iz_shifted-1)*x.PN[2] + iy_shifted-1)*x.PN[1] + ix_shifted

    if haskey(send_ranks,myrank_send)
    else
        send_ranks[myrank_send] = Data_sent_fermion(N,NC)
    end
    send_ranks[myrank_send].count += 1
    send_ranks[myrank_send].data[:,:,send_ranks[myrank_send].count] .= tempmatrix_mini
    send_ranks[myrank_send].positions[send_ranks[myrank_send].count] = disp

end

function mpi_updates_fermion_1data!(x::WilsonFermion_4D_nowing_mpi{NC},send_ranks) where NC
    if length(send_ranks) != 0
        NG = x.NG
        #=
        for rank=0:get_nprocs(U)
            if rank == get_myrank(U)
                println("myrank = ",myrank)
                for (key,value) in send_ranks
                    println(key,"\t",value.count)
                end
            end
            barrier(U)
        end
        =#
        tempmatrix = x.tempmatrix #zeros(ComplexF64,NC,NC,N)
        #tempmatrix = zeros(ComplexF64,NC,NC,N)
        positions = x.positions

        win = x.win
        #@time win = MPI.Win_create(tempmatrix,comm)
        #println(typeof(win))
        #Isend Irecv
        MPI.Win_fence(0, win)

        for (myrank_send,value) in send_ranks
            count = value.count
            MPI.Put(value.data[:,:,1:count], myrank_send,win)
        end

        MPI.Win_fence(0, win)
        #MPI.free(win)

        win_i = x.win_i#MPI.Win_create(positions,comm)
        MPI.Win_fence(0, win_i)

        for (myrank_send,value) in send_ranks
            count = value.count
            MPI.Put(value.positions[1:count], myrank_send,win_i)
        end

        MPI.Win_fence(0, win_i)
        #MPI.free(win_i)

        countvec = x.countvec#zeros(Int64,1)
        win_c = x.win_1i
        #win_c = MPI.Win_create(countvec,comm)
        MPI.Win_fence(0, win_c)

        for (myrank_send,value) in send_ranks
            count = value.count
            MPI.Put(Int64[count], myrank_send,win_c)
        end

        MPI.Win_fence(0, win_c)
        #MPI.free(win_c)

        count = countvec[1]

        #=
        for rank=0:get_nprocs(U)
            if rank == get_myrank(U)
                println("myrank = ",myrank)
                for position in positions[1:count]
                    println(position)
                end
            end
            barrier(U)
        end
        =#

        for i = 1:count
            position = positions[i]
            for jc = 1:NG
                for ic= 1:NC
                    ii = ((position-1)*NG+ jc-1)*NC + ic
                    x.fshifted[ii] = tempmatrix[ic,jc,i]
                end
            end
            #println(position)
        end

        #error("in shiftdU")
    end
end

function mpi_updates_fermion_moredata!(x::WilsonFermion_4D_nowing_mpi{NC},send_ranks) where {NC}
    NG = x.NG
    otherranks = x.otherranks
    win_other = x.win_other

    MPI.Win_fence(0, win_other)
    myrank = get_myrank(x)
    nprocs = get_nprocs(x)
    for (myrank_send,value) in send_ranks
        count = value.count
        MPI.Put(Int64[count], myrank_send,myrank,win_other)
    end
    MPI.Win_fence(0, win_other)

    tempmatrix = x.tempmatrix #zeros(ComplexF64,NC,NC,N)
    #tempmatrix = zeros(ComplexF64,NC,NC,N)
    positions = x.positions

    win = x.win
    #@time win = MPI.Win_create(tempmatrix,comm)
    #println(typeof(win))
    #Isend Irecv

    win_i = x.win_i#MPI.Win_create(positions,comm)

    win_c = x.win_1i
    #win_c = MPI.Win_create(countvec,comm)

    countvec = x.countvec#zeros(Int64,1)

    your_ranks = x.your_ranks #zeros(Int64,nprocs,nprocs)
    your_ranks .= -1

    MPI.Win_fence(0, win_other)
    icount = 0
    for (myrank_send,value) in send_ranks
        icount += 1
        MPI.Get(view(your_ranks,1:nprocs,icount), myrank_send,win_other)
    end
    MPI.Win_fence(0, win_other)

    MPI.Win_fence(0, win)
    MPI.Win_fence(0, win_i)
    MPI.Win_fence(0, win_c)

    icount = 0
    for (myrank_send,value) in send_ranks
        count = value.count
        icount += 1
        disp = 0
        for irank=1:myrank
            if your_ranks[irank,icount] != -1
                disp += your_ranks[irank,icount]
            end
        end

        MPI.Put(value.positions[1:count], myrank_send,disp,win_i)
        MPI.Put(value.data[:,:,1:count], myrank_send,disp*NC*NG,win)
    end

    MPI.Win_fence(0, win)
    MPI.Win_fence(0, win_i)
    MPI.Win_fence(0, win_c)

    your_ranks .= -1

    totaldatanum = sum(otherranks)

    for i = 1:totaldatanum
        position = positions[i]
        for jc = 1:NG
            for ic= 1:NC
                ii = ((position-1)*NG + jc-1)*NC + ic
                x.fshifted[ii] = tempmatrix[ic,jc,i]
            end
        end
        #println(position)
    end

    otherranks .= 0

end

function mpi_updates_fermion!(x::WilsonFermion_4D_nowing_mpi{NC},send_ranks) where {NC}
    if length(send_ranks) != 0

        val = MPI.Allreduce(length(send_ranks), +,x.comm) ÷ get_nprocs(x)

        #=
        for rank=0:get_nprocs(x)
            if rank == get_myrank(x)
                println("length = ",val,"\t")
                println("myrank = ",rank," length = $(length(send_ranks))")
            end
            barrier(x)
        end
        =#

        #if val == 1
        #    mpi_updates_fermion_1data!(x,send_ranks)
        #else
            mpi_updates_fermion_moredata!(x,send_ranks)
        #end

        return
    end
end

function shifted_fermion!(x::WilsonFermion_4D_nowing_mpi{NC},boundarycondition,shift) where {NC}
    PEs = x.PEs
    PN = x.PN
    myrank = x.myrank
    myrank_xyzt = x.myrank_xyzt
    myrank_xyzt_send = x.myrank_xyzt
    bc = boundarycondition
    NG = x.NG
    #tempmatrix = zeros(ComplexF64,NC,NC)#view(x.tempmatrix,1:NC,1:NC,1) #zeros(ComplexF64,NC,NC)
    tempmatrix_mini = view(x.tempmatrix,1:NC,1:NG,1)
    lat_size = size(x.fshifted)
    send_ranks = x.send_ranks
    empty!(send_ranks)
    # Dict{Int64,Data_sent}()
    N =  prod(x.PN)

    for it=1:x.PN[4]
        it_shifted = it - shift[4]
        it_global = myrank_xyzt[4]*x.PN[4] + it
        it_shifted_global = it_global - shift[4]
        inside_up = it_shifted_global >x. NT
        inside_down = it_shifted_global < 1
        factor_t = ifelse(inside_up || inside_down,bc[4],1)
        #if myrank_xyzt[4] == 0
        while it_shifted_global < 1
            it_shifted += x.NT
            it_shifted_global += x.NT
        end
            #it_shifted += ifelse(it_shifted < 1,x.NT,0)
        #end
        #if myrank_xyzt[4] == PEs[4]-1
        while it_shifted_global > x.NT
            it_shifted += -x.NT
            it_shifted_global += -x.NT
        end
            #it_shifted += ifelse(it_shifted > x.PN[4],-x.NT,0)
        #end
        if it_shifted <= 0
            tP = div(it_shifted,x.PN[4]) -1
        else
            tP = div(it_shifted-1,x.PN[4])
        end
        #if tP < 0
        #    println("it_shifted $it_shifted tP = $tP myrank_xyzt $myrank_xyzt it = $it shift = $shift it_shifted_global $it_shifted_global")
        #end

        #it_shifted += ifelse(it_shifted < 1,x.PN[4],0)
        while it_shifted < 1
            it_shifted += x.PN[4]
        end
        while it_shifted > x.PN[4]
            it_shifted += -x.PN[4]
        end
        #it_shifted += ifelse(it_shifted > x.PN[4],-x.PN[4],0)

        for iz=1:x.PN[3]
            iz_shifted = iz - shift[3]
            iz_global = myrank_xyzt[3]*x.PN[3] + iz
            iz_shifted_global = iz_global - shift[3]
            inside_up = iz_shifted_global >x. NZ
            inside_down = iz_shifted_global < 1
            factor_z = ifelse(inside_up || inside_down,bc[3],1)
            #if myrank_xyzt[3] == 0
            while iz_shifted_global < 1
                iz_shifted += x.NZ
                iz_shifted_global += x.NZ
            end
                #iz_shifted += ifelse(iz_shifted < 1,x.NZ,0)
            #end
            #if myrank_xyzt[3] == PEs[3]-1
            while iz_shifted_global  > x.NZ
                iz_shifted += -x.NZ
                iz_shifted_global += -x.NZ
            end

                #iz_shifted += ifelse(iz_shifted > x.PN[3],-x.NZ,0)
            #end

            if iz_shifted <=  0
                zP = div(iz_shifted,x.PN[3])-1
            else
                zP = div(iz_shifted-1,x.PN[3])
            end

            while iz_shifted < 1
                iz_shifted += x.PN[3]
            end
            while iz_shifted > x.PN[3]
                iz_shifted += -x.PN[3]
            end
            #iz_shifted += ifelse(iz_shifted < 1,x.PN[3],0)
            #iz_shifted += ifelse(iz_shifted > x.PN[3],-x.PN[3],0)

            for iy=1:x.PN[2]

                iy_shifted = iy - shift[2]
                iy_global = myrank_xyzt[2]*x.PN[2] + iy
                iy_shifted_global = iy_global - shift[2]
                inside_up = iy_shifted_global >x. NY
                inside_down = iy_shifted_global < 1
                factor_y = ifelse(inside_up || inside_down,bc[2],1)
                #if myrank_xyzt[2] == 0
                while iy_shifted_global < 1
                    iy_shifted += x.NY
                    iy_shifted_global += x.NY
                end

                    #iy_shifted += ifelse(iy_shifted < 1,x.NY,0)
                #end
                #if myrank_xyzt[2] == PEs[2]-1
                while iy_shifted_global > x.NY
                    iy_shifted += -x.NY
                    iy_shifted_global += -x.NY
                end
                    #iy_shifted += ifelse(iy_shifted > x.PN[2],-x.NY,0)
                #end

                if iy_shifted <=  0
                    yP = div(iy_shifted,x.PN[2])-1
                else
                    yP = div(iy_shifted-1,x.PN[2])
                end

                while iy_shifted < 1
                    iy_shifted += x.PN[2]
                end
                while iy_shifted > x.PN[2]
                    iy_shifted += -x.PN[2]
                end
                #iy_shifted += ifelse(iy_shifted < 1,x.PN[2],0)
                #iy_shifted += ifelse(iy_shifted > x.PN[2],-x.PN[2],0)

                for ix=1:x.PN[1]
                    ix_shifted = ix - shift[1]
                    ix_global = myrank_xyzt[1]*x.PN[1] + ix
                    ix_shifted_global = ix_global - shift[1]
                    inside_up = ix_shifted_global >x. NX
                    inside_down = ix_shifted_global < 1
                    factor_x = ifelse(inside_up || inside_down,bc[1],1)
                    #if myrank_xyzt[1] == 0
                    while ix_shifted_global < 1
                        ix_shifted += x.NX
                        ix_shifted_global += x.NX
                    end
                        #ix_shifted += ifelse(ix_shifted < 1,x.NX,0)
                    #end
                    #if myrank_xyzt[1] == PEs[1]-1
                    while ix_shifted_global > x.NX
                        ix_shifted += -x.NX
                        ix_shifted_global += -x.NX
                    end
                        #ix_shifted += ifelse(ix_shifted > x.PN[1],-x.NX,0)
                    #end

                    if ix_shifted <=  0
                        xP = div(ix_shifted,x.PN[1])-1
                    else
                        xP = div(ix_shifted-1,x.PN[1])
                    end

                    while ix_shifted < 1
                        ix_shifted += x.PN[1]
                    end
                    while ix_shifted > x.PN[1]
                        ix_shifted += -x.PN[1]
                    end
                    #ix_shifted += ifelse(ix_shifted < 1,x.PN[1],0)
                    #ix_shifted += ifelse(ix_shifted > x.PN[1],-x.PN[1],0)
                    #xP = div(ix_shifted-1,x.PN[1])
                    #println((tP,zP,yP,xP),"\t $shift")
                    if tP == 0 && zP == 0 && yP == 0 && xP == 0
                        for jc=1:NG
                            @simd for ic=1:NC
                                #v = getvalue(U,ic,jc,ix_shifted,iy_shifted,iz_shifted,it_shifted)
                                #x.Ushifted[ic,jc,ix,iy,iz,it] = v
                                v = factor_x*factor_y*factor_z*factor_t*getvalue(x,ic,jc,ix,iy,iz,it)
                                x.fshifted[ic,jc,ix_shifted,iy_shifted,iz_shifted,it_shifted] = v

                            end
                        end
                    else
                        update_sent_data!(send_ranks,N,ix,iy,iz,it,
                        ix_shifted,iy_shifted,iz_shifted,it_shifted,PEs,myrank_xyzt,xP,yP,zP,tP,x,
                        factor_x*factor_y*factor_z*factor_t)

                    end
                end
            end
        end
    end

    barrier(x)

    if length(send_ranks) != 0
        mpi_updates_fermion!(x,send_ranks)
    end

end

function setvalue!(F::T,v,i1,i2,i3,i4,i5,i6)  where T <: Shifted_fermionfields_4D_nowing_mpi
    error("type $(typeof(F)) has no setindex method. This type is read only.")
end

function getvalue(F::T,i1,i2,i3,i4,i5,i6)  where T <: Shifted_fermionfields_4D_nowing
    @inbounds return F.parent.fshifted[i1,i2 ,i3 ,i4 ,i5,i6]
end

"""
c-------------------------------------------------c
c     Random number function for Gaussian  Noise
    with σ^2 = 1/2
c-------------------------------------------------c
    """
function gauss_distribution_fermion!(x::WilsonFermion_4D_nowing_mpi{NC}) where {NC}
    NX = x.NX
    NY = x.NY
    NZ = x.NZ
    NT = x.NT
    NG = x.NG
    #n6 = size(x.f)[6]
    σ = sqrt(1/2)

    for it=1:x.PN[4]
        for iz=1:x.PN[3]
            for iy=1:x.PN[2]
                for ix=1:x.PN[1]
                    for ialpha = 1:NG
                        @inbounds @simd for ic=1:NC
                            v = σ*randn()+im*σ*randn()
                            setvalue!(x,v,ic,ialpha,ix,iy,iz,it)
                            #x[ic,ialpha,ix,iy,iz,it] = σ*randn()+im*σ*randn()
                        end
                    end
                end
            end
        end
    end
    set_wing_fermion!(x)
    return
end

"""
c-------------------------------------------------c
c     Random number function for Gaussian  Noise
    with σ^2 = 1/2
c-------------------------------------------------c
    """
function gauss_distribution_fermion!(x::WilsonFermion_4D_nowing_mpi{NC},randomfunc,σ) where {NC}

    NX = x.NX
    NY = x.NY
    NZ = x.NZ
    NT = x.NT
    NG = x.NG
    #n6 = size(x.f)[6]
    #σ = sqrt(1/2)

    for it=1:x.PN[4]
        for iz=1:x.PN[3]
            for iy=1:x.PN[2]
                for ix=1:x.PN[1]
                    for mu = 1:NG
                        @inbounds @simd for ic=1:NC

                            v1 = sqrt(-log(randomfunc()+1e-10))
                            v2 = 2pi*randomfunc()

                            xr = v1*cos(v2)
                            xi = v1 * sin(v2)

                            v = σ*xr + σ*im*xi

                            setvalue!(x,v,ic,mu,ix,iy,iz,it)

                            #x[ic,ix,iy,iz,it,mu] = σ*xr + σ*im*xi
                        end
                    end
                end
            end
        end
    end
    #error("v")

    set_wing_fermion!(x)

    return
end

function Z2_distribution_fermion!(x::AbstractFermionfields_4D{NC})  where NC
    NX = x.NX
    NY = x.NY
    NZ = x.NZ
    NT = x.NT
    #n6 = size(x.f)[6]
    #σ = sqrt(1/2)

    for it=1:x.PN[4]
        for iz=1:x.PN[3]
            for iy=1:x.PN[2]
                for ix=1:x.PN[1]
                    for mu = 1:NG
                        for ic=1:NC
                            v = rand([-1,1])
                            setvalue!(x,v,ic,mu,ix,iy,iz,it)
                            #x[ic,ix,iy,iz,it,mu] = rand([-1,1])
                        end
                    end
                end
            end
        end
    end

    set_wing_fermion!(x)

    return
end

"""
c-------------------------------------------------c
c     Random number function Z4  Noise
c     https://arxiv.org/pdf/1611.01193.pdf
c-------------------------------------------------c
    """
function Z4_distribution_fermi!(x::AbstractFermionfields_4D{NC})  where NC
        NX = x.NX
        NY = x.NY
        NZ = x.NZ
        NT = x.NT
        n6 = size(x.f)[6]
        θ = 0.0
        N::Int32 = 4
        Ninv = Float64(1/N)
        for it=1:x.PN[4]
            for iz=1:x.PN[3]
                for iy=1:x.PN[2]
                    for ix=1:x.PN[1]
                        for ialpha = 1:NG
                            @inbounds @simd for ic=1:NC
                                θ = Float64(rand(0:N-1))*π*Ninv # r \in [0,π/4,2π/4,3π/4]
                                v = cos(θ)+im*sin(θ)
                                setvalue!(x,v,ic,ialpha,ix,iy,iz,it)
                                #x[ic,ix,iy,iz,it,ialpha] = cos(θ)+im*sin(θ)
                            end
                        end
                    end
                end
            end
        end

        set_wing_fermion!(x)

        return
end

function gauss_distribution_fermion!(x::AbstractFermionfields_4D{NC},randomfunc) where NC
    σ = 1
    gauss_distribution_fermion!(x,randomfunc,σ)
end

function set_wing_fermion!(a::WilsonFermion_4D_nowing_mpi{NC},boundarycondition) where {NC}
    return
end

function Wx!(xout::T,U::Array{G,1},x::T,A)  where  {T <: WilsonFermion_4D_nowing_mpi,G <: AbstractGaugefields}
    #temps::Array{T,1},boundarycondition) where  {T <: WilsonFermion_4D_nowing_mpi,G <: AbstractGaugefields}
    temp = A._temporary_fermi[4]#temps[4]
    temp1 = A._temporary_fermi[1] #temps[1]
    temp2 = A._temporary_fermi[2] #temps[2]

    #temp = temps[4]
    #temp1 = temps[1]
    #temp2 = temps[2]

    clear_fermion!(temp)
    #set_wing_fermion!(x)
    for ν=1:4

        xplus = shift_fermion(x,ν)
        #println(xplus)

        mul!(temp1,U[ν],xplus)

        #fermion_shift!(temp1,U,ν,x)

        #... Dirac multiplication

        mul!(temp1,view(A.rminusγ,:,:,ν))

        xminus = shift_fermion(x,-ν)
        Uminus = shift_U(U[ν],-ν)

        mul!(temp2,Uminus',xminus)

        #
        #fermion_shift!(temp2,U,-ν,x)
        #mul!(temp2,view(x.rplusγ,:,:,ν),temp2)
        mul!(temp2,view(A.rplusγ,:,:,ν))

        add_fermion!(temp,A.hopp[ν],temp1,A.hopm[ν],temp2)

    end

    clear_fermion!(xout)
    add_fermion!(xout,1,x,-1,temp)

    set_wing_fermion!(xout,A.boundarycondition)

    #display(xout)
    #    exit()
    return
end

function Wdagx!(xout::T,U::Array{G,1},
    x::T,A) where  {T <: WilsonFermion_4D_nowing_mpi,G <: AbstractGaugefields}
    #,temps::Array{T,1},boundarycondition) where  {T <: WilsonFermion_4D_nowing_mpi,G <: AbstractGaugefields}
    temp = A._temporary_fermi[4] #temps[4]
    temp1 = A._temporary_fermi[1] #temps[1]
    temp2 = A._temporary_fermi[2] #temps[2]

    clear_fermion!(temp)
    #set_wing_fermion!(x)
    for ν=1:4
        xplus = shift_fermion(x,ν)
        mul!(temp1,U[ν],xplus)

        #fermion_shift!(temp1,U,ν,x)

        #... Dirac multiplication
        #mul!(temp1,view(x.rminusγ,:,:,ν),temp1)
        mul!(temp1,view(A.rplusγ,:,:,ν))

        #
        xminus = shift_fermion(x,-ν)
        Uminus = shift_U(U[ν],-ν)

        mul!(temp2,Uminus',xminus)
        #fermion_shift!(temp2,U,-ν,x)
        #mul!(temp2,view(x.rminusγ,:,:,ν),temp2)
        mul!(temp2,view(A.rminusγ,:,:,ν))

        add_fermion!(temp,A.hopp[ν],temp1,A.hopm[ν],temp2)

    end

    clear_fermion!(xout)
    add_fermion!(xout,1,x,-1,temp)
    set_wing_fermion!(xout,A.boundarycondition)

    #display(xout)
    #    exit()
    return
end

function add_fermion!(c::WilsonFermion_4D_nowing_mpi{NC},α::Number,a::T1,β::Number,b::T2) where {NC,T1 <: Abstractfermion,T2 <: Abstractfermion}#c += alpha*a + beta*b
    #n1,n2,n3,n4,n5,n6 = size(c.f)

    @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 ialpha=1:4
                        @simd for k1=1:NC
                            v = getvalue(c,k1,ialpha,ix,iy,iz,it) + α*getvalue(a,k1,ialpha,ix,iy,iz,it)+β*getvalue(b,k1,ialpha,ix,iy,iz,it)
                                setvalue!(c,v,k1,ialpha,ix,iy,iz,it)
                            #println(a.f[i1,i2,i3,i4,i5,i6],"\t",b.f[i1,i2,i3,i4,i5,i6] )
                            #c.f[i1,i2,i3,i4,i5,i6] += α*a.f[i1,i2,i3,i4,i5,i6] + β*b.f[i1,i2,i3,i4,i5,i6]
                        end
                    end
                end
            end
        end
    end
    return
end

function add_fermion!(c::WilsonFermion_4D_nowing_mpi{NC},α::Number,a::T1) where {NC,T1 <: Abstractfermion}#c += alpha*a
    #n1,n2,n3,n4,n5,n6 = size(c.f)

    @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 ialpha=1:4
                        @simd for k1=1:NC
                            v = getvalue(c,k1,ialpha,ix,iy,iz,it) + α*getvalue(a,k1,ialpha,ix,iy,iz,it)
                            setvalue!(c,v,k1,ialpha,ix,iy,iz,it)
                            #println(a.f[i1,i2,i3,i4,i5,i6],"\t",b.f[i1,i2,i3,i4,i5,i6] )
                            #c.f[i1,i2,i3,i4,i5,i6] += α*a.f[i1,i2,i3,i4,i5,i6]
                        end
                    end
                end
            end
        end
    end
    return
end

function add_fermion!(c::WilsonFermion_4D_nowing_mpi{NC},α::Number,a::T1,β::Number,b::T2,iseven) where {NC,T1 <: Abstractfermion,T2 <: Abstractfermion}#c += alpha*a + beta*b

    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
                    evenodd = ifelse((ix + iy + iz + it) % 2 == 0,true,false)
                    if evenodd == iseven
                        for ialpha=1:4
                            @simd for k1=1:NC
                                v = getvalue(c,k1,ialpha,ix,iy,iz,it) + α*getvalue(a,k1,ialpha,ix,iy,iz,it)+β*getvalue(b,k1,ialpha,ix,iy,iz,it)
                                setvalue!(c,v,k1,ialpha,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_fermion!(c,iseven)

    return
end

function add_fermion!(c::WilsonFermion_4D_nowing_mpi{NC},α::Number,a::T1,iseven::Bool) where {NC,T1 <: Abstractfermion,T2 <: Abstractfermion}#c += alpha*a + beta*b

    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
                        evenodd = ifelse((ix + iy + iz + it) % 2 == 0,true,false)
                        if evenodd == iseven
                            for ialpha=1:4
                            @simd for k1=1:NC
                                v = getvalue(c,k1,ialpha,ix,iy,iz,it) + α*getvalue(a,k1,ialpha,ix,iy,iz,it)
                                setvalue!(c,v,k1,ialpha,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_fermion!(c,iseven)

    return
end

function WWx!(xout::T,U::Array{G,1},x::T,A)  where  {T <: WilsonFermion_4D_nowing_mpi,G <: AbstractGaugefields} #(1 - K^2 Teo Toe) xe
    iseven = true
    temp = A._temporary_fermi[4]#temps[4]
    temp2 = A._temporary_fermi[5]#temps[4]
    Toex!(temp2,U,x,A,iseven)

    iseven = false
    Toex!(temp,U,temp2,A,iseven)
    add_fermion!(xout,1,x,-1,temp)

    iseven = true
    set_wing_fermion!(xout,A.boundarycondition)

    return
end

function clear_fermion!(a::WilsonFermion_4D_nowing_mpi{NC} ,iseven) where {NC}
    n1,n6,n2,n3,n4,n5= size(a.f)
    @inbounds for i5=1:n5
        it = i5
        for i4=1:n4
            iz = i4
            for i3=1:n3
                iy = i3
                for i2=1:n2
                    ix = i2
                    evenodd = ifelse((ix+iy+iz+it) % 2 == 0,true,false)
                    if evenodd == iseven
                        for i6=1:n6
                            @simd for i1=1:NC
                                a.f[i1,i6,i2,i3,i4,i5]= 0
                            end
                        end
                    end
                end
            end
        end
    end
end

function LinearAlgebra.dot(a::WilsonFermion_4D_nowing_mpi{NC},b::WilsonFermion_4D_nowing_mpi{NC}) where {NC}
    NT = a.NT
    NZ = a.NZ
    NY = a.NY
    NX = a.NX
    NG = a.NG

    c = 0.0im
    @inbounds  for it=1:a.PN[4]
        for iz=1:a.PN[3]
            for iy=1:a.PN[2]
                for ix=1:a.PN[1]
                    for α=1:NG
                        @simd for ic=1:NC
                            va = getvalue(a,ic,α,ix,iy,iz,it)
                            vb = getvalue(b,ic,α,ix,iy,iz,it)
                            c += conj(va)*vb
                            #c+= conj(a[ic,ix,iy,iz,it,α])*b[ic,ix,iy,iz,it,α]
                        end
                    end
                end
            end
        end
    end

    c = MPI.Allreduce(c,MPI.SUM,a.comm)
    return c
end

function LinearAlgebra.mul!(y::WilsonFermion_4D_nowing_mpi{3},A::T,x::T3) where {T<:Abstractfields,T3 <:Abstractfermion}
    #@assert 3 == x.NC "dimension mismatch! NC in y is 3 but NC in x is $(x.NC)"
    NX = y.NX
    NY = y.NY
    NZ = y.NZ
    NT = y.NT
    NG = y.NG

    @inbounds for it=1:y.PN[4]
        for iz=1:y.PN[3]
            for iy=1:y.PN[2]
                for ix=1:y.PN[1]
                    for ialpha=1:NG
                        #println(ix)
                        x1 = getvalue(x,1,ialpha,ix,iy,iz,it)#  x[ic,1,ix,iy,iz,it)
                        x2 = getvalue(x,2,ialpha,ix,iy,iz,it)
                        x3 = getvalue(x,3,ialpha,ix,iy,iz,it)

                        v= getvalue(A,1,1,ix,iy,iz,it)*x1 +
                                                    getvalue(A,1,2,ix,iy,iz,it)*x2+
                                                    getvalue(A,1,3,ix,iy,iz,it)*x3
                        setvalue!(y,v,1,ialpha,ix,iy,iz,it)
                        v = getvalue(A,2,1,ix,iy,iz,it)*x1+
                                                    getvalue(A,2,2,ix,iy,iz,it)*x2 +
                                                    getvalue(A,2,3,ix,iy,iz,it)*x3
                        setvalue!(y,v,2,ialpha,ix,iy,iz,it)
                        v= getvalue(A,3,1,ix,iy,iz,it)*x1+
                                                    getvalue(A,3,2,ix,iy,iz,it)*x2 +
                                                    getvalue(A,3,3,ix,iy,iz,it)*x3
                        setvalue!(y,v,3,ialpha,ix,iy,iz,it)
                        # =#
                    end
                end
            end
        end
    end
end

function LinearAlgebra.mul!(y::WilsonFermion_4D_nowing_mpi{2},A::T,x::T3) where {T<:Abstractfields,T3 <:Abstractfermion}
    #@assert 3 == x.NC "dimension mismatch! NC in y is 3 but NC in x is $(x.NC)"
    NX = y.NX
    NY = y.NY
    NZ = y.NZ
    NT = y.NT
    NG = y.NG

    @inbounds for it=1:y.PN[4]
        for iz=1:y.PN[3]
            for iy=1:y.PN[2]
                for ix=1:y.PN[1]
                    for ialpha=1:NG
                        #println(ix)
                        x1 = getvalue(x,1,ialpha,ix,iy,iz,it)#  x[ic,1,ix,iy,iz,it)
                        x2 = getvalue(x,2,ialpha,ix,iy,iz,it)
                        #x3 = getvalue(x,3,ialpha,ix,iy,iz,it)

                        v= getvalue(A,1,1,ix,iy,iz,it)*x1 +
                                                    getvalue(A,1,2,ix,iy,iz,it)*x2#+
                                                    #getvalue(A,1,3,ix,iy,iz,it)*x3
                        setvalue!(y,v,1,ialpha,ix,iy,iz,it)
                        v = getvalue(A,2,1,ix,iy,iz,it)*x1+
                                                    getvalue(A,2,2,ix,iy,iz,it)*x2 #+
                                                    #getvalue(A,2,3,ix,iy,iz,it)*x3
                        setvalue!(y,v,2,ialpha,ix,iy,iz,it)
                        # =#
                    end
                end
            end
        end
    end
end

function LinearAlgebra.mul!(x::WilsonFermion_4D_nowing_mpi{NC},A::TA) where {TA <: AbstractMatrix, NC}
    NX = x.NX
    NY = x.NY
    NZ = x.NZ
    NT = x.NT

    #n6 = size(x.f)[6]
    #f = zeros(ComplexF64,4)
    #e = zeros(ComplexF64,4)

    @inbounds  for it=1:x.PN[4]
        for iz=1:x.PN[3]
            for iy=1:x.PN[2]
                for ix=1:x.PN[1]
                    for ic=1:NC
                            e1 = getvalue(x,ic,1,ix,iy,iz,it)#  x[ic,1,ix,iy,iz,it)
                            e2 = getvalue(x,ic,2,ix,iy,iz,it)
                            e3 = getvalue(x,ic,3,ix,iy,iz,it)
                            e4 = getvalue(x,ic,4,ix,iy,iz,it)

                            v = A[1,1]*e1+A[1,2]*e2+A[1,3]*e3+A[1,4]*e4
                            setvalue!(x,v,ic,1,ix,iy,iz,it)
                            v = A[2,1]*e1+A[2,2]*e2+A[2,3]*e3+A[2,4]*e4
                            setvalue!(x,v,ic,2,ix,iy,iz,it)
                            v = A[3,1]*e1+A[3,2]*e2+A[3,3]*e3+A[3,4]*e4
                            setvalue!(x,v,ic,3,ix,iy,iz,it)
                            v = A[4,1]*e1+A[4,2]*e2+A[4,3]*e3+A[4,4]*e4
                            setvalue!(x,v,ic,4,ix,iy,iz,it)

                            #x[ic,1,ix,iy,iz,it) = A[1,1]*e1+A[1,2]*e2+A[1,3]*e3+A[1,4]*e4
                            #x[ic,2,ix,iy,iz,it] = A[2,1]*e1+A[2,2]*e2+A[2,3]*e3+A[2,4]*e4
                            #x[ic,3,ix,iy,iz,it) = A[3,1]*e1+A[3,2]*e2+A[3,3]*e3+A[3,4]*e4
                            #x[ic,4,ix,iy,iz,it) = A[4,1]*e1+A[4,2]*e2+A[4,3]*e3+A[4,4]*e4

                    end
                end
            end
        end
    end

end

function LinearAlgebra.mul!(u::T1,x::Abstractfermion,y::Adjoint_fermionfields{<: WilsonFermion_4D_nowing_mpi{NC}}) where {T1 <: AbstractGaugefields,NC}
    #_,NX,NY,NZ,NT,NG = size(y)
    NG = x.NG
    clear_U!(u)

    @inbounds for it=1:x.PN[4]
        for iz=1:x.PN[3]
            for iy=1:x.PN[2]
                for ix=1:x.PN[1]
                    for ib=1:NC
                        for ik=1:NG
                            @simd for ia=1:NC
                                v = getvalue(u,ia,ib,ix,iy,iz,it) + getvalue(x,ia,ik,ix,iy,iz,it)*getvalue(y,ib,ik,ix,iy,iz,it)
                                setvalue!(u,v,ia,ib,ix,iy,iz,it)
                                #u[ia,ib,ix,iy,iz,it] += x[ia,ix,iy,iz,it,ik]*y[ib,ix,iy,iz,it,ik]
                            end
                        end
                    end
                end
            end
        end
    end
    set_wing_U!(u)
end

"""
mul!(u,x,y) -> u_{ab} = x_a*y_b
"""
function LinearAlgebra.mul!(u::T1,x::WilsonFermion_4D_nowing_mpi{NC},y::WilsonFermion_4D_nowing_mpi{NC}) where {T1 <: AbstractGaugefields,NC}
    NX = x.NX
    NY = x.NY
    NZ = x.NZ
    NT = x.NT
    NG = x.NG
    clear_U!(u)

    @inbounds  for it=1:x.PN[4]
        for iz=1:x.PN[3]
            for iy=1:x.PN[2]
                for ix=1:x.PN[1]
                    for ik=1:NG
                        for ib=1:NC
                            @simd for ia=1:NC
                                v = getvalue(u,ia,ib,ix,iy,iz,it) + getvalue(x,ia,ik,ix,iy,iz,it)*getvalue(y,ib,ik,ix,iy,iz,it)
                                setvalue!(u,v,ia,ib,ix,iy,iz,it)

                                #u[ia,ib,ix,iy,iz,it] += x[ia,ix,iy,iz,it,ik]*y[ib,ix,iy,iz,it,ik]
                            end
                        end
                    end
                end
            end
        end
    end
    set_wing_U!(u)
end

function LinearAlgebra.mul!(y::WilsonFermion_4D_nowing_mpi{3},x::T3,A::T) where {T<:Abstractfields,T3 <:Abstractfermion}
    #@assert 3 == x.NC "dimension mismatch! NC in y is 3 but NC in x is $(x.NC)"
    NX = y.NX
    NY = y.NY
    NZ = y.NZ
    NT = y.NT
    NG = y.NG

    @inbounds for it=1:y.PN[4]
        for iz=1:y.PN[3]
            for iy=1:y.PN[2]
                for ix=1:y.PN[1]
                    for ialpha=1:NG
                        x1 = getvalue(x,1,ialpha,ix,iy,iz,it)
                        x2 = getvalue(x,2,ialpha,ix,iy,iz,it)
                        x3 = getvalue(x,3,ialpha,ix,iy,iz,it)
                        v = x1*getvalue(A,1,1,ix,iy,iz,it) +
                                                    x2*getvalue(A,2,1,ix,iy,iz,it)+
                                                    x3*getvalue(A,3,1,ix,iy,iz,it)
                        setvalue!(y,v,1,ialpha,ix,iy,iz,it)
                        v = x1*getvalue(A,1,2,ix,iy,iz,it)+
                                                    x2*getvalue(A,2,2,ix,iy,iz,it) +
                                                    x3*getvalue(A,3,2,ix,iy,iz,it)
                        setvalue!(y,v,2,ialpha,ix,iy,iz,it)
                        v = x1*getvalue(A,1,3,ix,iy,iz,it)+
                                                    x2*getvalue(A,2,3,ix,iy,iz,it) +
                                                    x3*getvalue(A,3,3,ix,iy,iz,it)
                        setvalue!(y,v,3,ialpha,ix,iy,iz,it)
                    end
                end
            end
        end
    end
end

function LinearAlgebra.mul!(y::WilsonFermion_4D_nowing_mpi{NC},A::T,x::T3) where {NC,T<:Number,T3 <:Abstractfermion}
    @assert NC == x.NC "dimension mismatch! NC in y is $NC but NC in x is $(x.NC)"
    NX = y.NX
    NY = y.NY
    NZ = y.NZ
    NT = y.NT
    NG = y.NG

    @inbounds for it=1:y.PN[4]
        for iz=1:y.PN[3]
            for iy=1:y.PN[2]
                for ix=1:y.PN[1]
                    for ialpha=1:NG
                        for k1=1:NC
                            v = A*getvalue(x,k1,ialpha,ix,iy,iz,it)
                            setvalue!(y,v,k1,ialpha,ix,iy,iz,it) # A*getvalue(x,k1,ialpha,ix,iy,iz,it)
                        end
                    end
                end
            end
        end
    end
end

function LinearAlgebra.mul!(x::WilsonFermion_4D_nowing_mpi{NC},A::TA,iseven::Bool) where {TA <: AbstractMatrix, NC}
    NX = x.NX
    NY = x.NY
    NZ = x.NZ
    NT = x.NT

    #n6 = size(x.f)[6]
    #f = zeros(ComplexF64,4)
    #e = zeros(ComplexF64,4)

    @inbounds  for it=1:x.PN[4]
        for iz=1:x.PN[3]
            for iy=1:x.PN[2]
                for ix=1:x.PN[1]
                    evenodd = ifelse((ix+iy+iz+it) % 2 == 0,true,false)
                    if evenodd == iseven
                        for ic=1:NC
                            e1 = getvalue(x,ic,1,ix,iy,iz,it)#  x[ic,1,ix,iy,iz,it)
                            e2 = getvalue(x,ic,2,ix,iy,iz,it)
                            e3 = getvalue(x,ic,3,ix,iy,iz,it)
                            e4 = getvalue(x,ic,4,ix,iy,iz,it)

                            v = A[1,1]*e1+A[1,2]*e2+A[1,3]*e3+A[1,4]*e4
                            setvalue!(x,v,ic,1,ix,iy,iz,it)
                            v = A[2,1]*e1+A[2,2]*e2+A[2,3]*e3+A[2,4]*e4
                            setvalue!(x,v,ic,2,ix,iy,iz,it)
                            v = A[3,1]*e1+A[3,2]*e2+A[3,3]*e3+A[3,4]*e4
                            setvalue!(x,v,ic,3,ix,iy,iz,it)
                            v = A[4,1]*e1+A[4,2]*e2+A[4,3]*e3+A[4,4]*e4
                            setvalue!(x,v,ic,4,ix,iy,iz,it)
                        end

                        #x[ic,1,ix,iy,iz,it) = A[1,1]*e1+A[1,2]*e2+A[1,3]*e3+A[1,4]*e4
                        #x[ic,2,ix,iy,iz,it] = A[2,1]*e1+A[2,2]*e2+A[2,3]*e3+A[2,4]*e4
                        #x[ic,3,ix,iy,iz,it) = A[3,1]*e1+A[3,2]*e2+A[3,3]*e3+A[3,4]*e4
                        #x[ic,4,ix,iy,iz,it) = A[4,1]*e1+A[4,2]*e2+A[4,3]*e3+A[4,4]*e4
                    end

                end
            end
        end
    end

end

function LinearAlgebra.mul!(xout::WilsonFermion_4D_nowing_mpi{NC},A::TA,x::WilsonFermion_4D_nowing_mpi{NC}) where {TA <: AbstractMatrix, NC}
    NX = x.NX
    NY = x.NY
    NZ = x.NZ
    NT = x.NT

    #n6 = size(x.f)[6]
    #f = zeros(ComplexF64,4)
    #e = zeros(ComplexF64,4)
    @inbounds  for it=1:x.PN[4]
        for iz=1:x.PN[3]
            for iy=1:x.PN[2]
                for ix=1:x.PN[1]
                    for ic=1:NC
                            e1 = getvalue(x,ic,1,ix,iy,iz,it)#  x[ic,ix,iy,iz,it,1]
                            e2 = getvalue(x,ic,2,ix,iy,iz,it)
                            e3 = getvalue(x,ic,3,ix,iy,iz,it)
                            e4 = getvalue(x,ic,4,ix,iy,iz,it)

                            v = A[1,1]*e1+A[1,2]*e2+A[1,3]*e3+A[1,4]*e4
                            setvalue!(xout,v,ic,1,ix,iy,iz,it)
                            v = A[2,1]*e1+A[2,2]*e2+A[2,3]*e3+A[2,4]*e4
                            setvalue!(xout,v,ic,2,ix,iy,iz,it)
                            v = A[3,1]*e1+A[3,2]*e2+A[3,3]*e3+A[3,4]*e4
                            setvalue!(xout,v,ic,3,ix,iy,iz,it)
                            v = A[4,1]*e1+A[4,2]*e2+A[4,3]*e3+A[4,4]*e4
                            setvalue!(xout,v,ic,4,ix,iy,iz,it)
                    end
                end
            end
        end
    end

end

function LinearAlgebra.mul!(xout::WilsonFermion_4D_nowing_mpi{NC},A::TA,x::WilsonFermion_4D_nowing_mpi{NC},iseven::Bool) where {TA <: AbstractMatrix, NC}
    NX = x.NX
    NY = x.NY
    NZ = x.NZ
    NT = x.NT

    #n6 = size(x.f)[6]
    #f = zeros(ComplexF64,4)
    #e = zeros(ComplexF64,4)

    @inbounds  for it=1:x.PN[4]
        for iz=1:x.PN[3]
            for iy=1:x.PN[2]
                for ix=1:x.PN[1]
                    evenodd = ifelse((ix+iy+iz+it) % 2 == 0,true,false)
                    if evenodd == iseven
                        for ic=1:NC

                            e1 = getvalue(x,ic,1,ix,iy,iz,it)#  x[ic,ix,iy,iz,it,1]
                            e2 = getvalue(x,ic,2,ix,iy,iz,it)
                            e3 = getvalue(x,ic,3,ix,iy,iz,it)
                            e4 = getvalue(x,ic,4,ix,iy,iz,it)

                            v = A[1,1]*e1+A[1,2]*e2+A[1,3]*e3+A[1,4]*e4
                            setvalue!(xout,v,ic,ix,iy,iz,it,1)
                            v = A[2,1]*e1+A[2,2]*e2+A[2,3]*e3+A[2,4]*e4
                            setvalue!(xout,v,ic,2,ix,iy,iz,it)
                            v = A[3,1]*e1+A[3,2]*e2+A[3,3]*e3+A[3,4]*e4
                            setvalue!(xout,v,ic,3,ix,iy,iz,it)
                            v = A[4,1]*e1+A[4,2]*e2+A[4,3]*e3+A[4,4]*e4
                            setvalue!(xout,v,ic,4,ix,iy,iz,it)
                        end
                    end

                end
            end
        end
    end

end

function LinearAlgebra.mul!(xout::WilsonFermion_4D_nowing_mpi{NC},x::WilsonFermion_4D_nowing_mpi{NC},A::TA) where {TA <: AbstractMatrix, NC}
    NX = x.NX
    NY = x.NY
    NZ = x.NZ
    NT = x.NT

    #n6 = size(x.f)[6]
    #f = zeros(ComplexF64,4)
    #e = zeros(ComplexF64,4)

    @inbounds  for it=1:x.PN[4]
        for iz=1:x.PN[3]
            for iy=1:x.PN[2]
                for ix=1:x.PN[1]
                    for ic=1:NC
                            e1 = getvalue(x,ic,1,ix,iy,iz,it)#  x[ic,1,ix,iy,iz,it)
                            e2 = getvalue(x,ic,2,ix,iy,iz,it)
                            e3 = getvalue(x,ic,3,ix,iy,iz,it)
                            e4 = getvalue(x,ic,4,ix,iy,iz,it)

                            v = A[1,1]*e1+A[2,1]*e2+A[3,1]*e3+A[4,1]*e4
                            setvalue!(xout,v,ic,1,ix,iy,iz,it)
                            v = A[1,2]*e1+A[2,2]*e2+A[3,2]*e3+A[4,2]*e4
                            setvalue!(xout,v,ic,2,ix,iy,iz,it)
                            v = A[1,3]*e1+A[2,3]*e2+A[3,3]*e3+A[4,3]*e4
                            setvalue!(xout,v,ic,3,ix,iy,iz,it)
                            v = A[1,4]*e1+A[2,4]*e2+A[3,4]*e3+A[4,4]*e4
                            setvalue!(xout,v,ic,4,ix,iy,iz,it)

                            #=
                            e1 = x[ic,1,ix,iy,iz,it)
                            e2 = x[ic,2,ix,iy,iz,it]
                            e3 = x[ic,3,ix,iy,iz,it)
                            e4 = x[ic,4,ix,iy,iz,it)

                            xout[ic,1,ix,iy,iz,it) = A[1,1]*e1+A[2,1]*e2+A[3,1]*e3+A[4,1]*e4
                            xout[ic,2,ix,iy,iz,it] = A[1,2]*e1+A[2,2]*e2+A[3,2]*e3+A[4,2]*e4
                            xout[ic,3,ix,iy,iz,it) = A[1,3]*e1+A[2,3]*e2+A[3,3]*e3+A[4,3]*e4
                            xout[ic,4,ix,iy,iz,it) = A[1,4]*e1+A[2,4]*e2+A[3,4]*e3+A[4,4]*e4
                            =#

                    end
                end
            end
        end
    end

end

function LinearAlgebra.mul!(xout::WilsonFermion_4D_nowing_mpi{NC},x::WilsonFermion_4D_nowing_mpi{NC},A::TA,iseven) where {TA <: AbstractMatrix, NC}
    NX = x.NX
    NY = x.NY
    NZ = x.NZ
    NT = x.NT

    #n6 = size(x.f)[6]
    #f = zeros(ComplexF64,4)
    #e = zeros(ComplexF64,4)

    @inbounds  for it=1:x.PN[4]
        for iz=1:x.PN[3]
            for iy=1:x.PN[2]
                for ix=1:x.PN[1]
                    evenodd = ifelse((ix+iy+iz+it) % 2 == 0,true,false)
                    if evenodd == iseven
                        for ic=1:NC
                            e1 = getvalue(x,ic,1,ix,iy,iz,it)#  x[ic,1,ix,iy,iz,it)
                            e2 = getvalue(x,ic,2,ix,iy,iz,it)
                            e3 = getvalue(x,ic,3,ix,iy,iz,it)
                            e4 = getvalue(x,ic,4,ix,iy,iz,it)

                            v = A[1,1]*e1+A[2,1]*e2+A[3,1]*e3+A[4,1]*e4
                            setvalue!(xout,v,ic,1,ix,iy,iz,it)
                            v = A[1,2]*e1+A[2,2]*e2+A[3,2]*e3+A[4,2]*e4
                            setvalue!(xout,v,ic,2,ix,iy,iz,it)
                            v = A[1,3]*e1+A[2,3]*e2+A[3,3]*e3+A[4,3]*e4
                            setvalue!(xout,v,ic,3,ix,iy,iz,it)
                            v = A[1,4]*e1+A[2,4]*e2+A[3,4]*e3+A[4,4]*e4
                            setvalue!(xout,v,ic,4,ix,iy,iz,it)

                        end
                    end
                end
            end
        end
    end

end

"""
c--------------------------------------------------------------------------c
c     y = gamma_5 * x
c     here
c                  ( -1       )
c        GAMMA5 =  (   -1     )
c                  (     +1   )
c                  (       +1 )
c--------------------------------------------------------------------------c
    """
    function mul_γ5x!(y::WilsonFermion_4D_nowing_mpi{NC},x::WilsonFermion_4D_nowing_mpi{NC}) where {NC}
        n1,n6,n2,n3,n4,n5 = size(x.f)
        #println("axpby")

        @inbounds for i5=1:n5

            for i4=1:n4

                for i3=1:n3

                    for i2=1:n2

                        for i6=1:n6
                            @simd for i1=1:NC
                                y.f[i1,i6,i2,i3,i4,i5]=x.f[i1,i6,i2,i3,i4,i5]*ifelse(i6 <= 2,-1,1)
                            end
                        end
                    end
                end
            end
        end

    end

    function apply_γ5!(x::WilsonFermion_4D_nowing_mpi{NC}) where {NC}
        n1,n6,n2,n3,n4,n5 = size(x.f)
        #println("axpby")

        @inbounds for i5=1:n5
            #it = i5+NDW
            for i4=1:n4
                #iz = i4+NDW
                for i3=1:n3
                    #iy = i3+NDW
                    for i2=1:n2
                        #ix = i2+NDW
                        for i6=1:n6
                            @simd for i1=1:NC
                                x.f[i1,i6,i2,i3,i4,i5]=x.f[i1,i6,i2,i3,i4,i5]*ifelse(i6 <= 2,-1,1)
                            end
                        end
                    end
                end
            end
        end

    end

    function mul_1plusγ5x!(y::WilsonFermion_4D_nowing_mpi{NC},x::WilsonFermion_4D_nowing_mpi{NC})  where {NC}#(1+gamma_5)/2
        n1,n6,n2,n3,n4,n5 = size(x.f)
        #println("axpby")

        #@inbounds for i6=1:n6
            @inbounds for i5=1:n5
                #it = i5+NDW
                for i4=1:n4
                    #iz = i4+NDW
                    for i3=1:n3
                        #iy = i3+NDW
                        for i2=1:n2
                            #ix = i2+NDW
                            @simd for i1=1:NC

                                y.f[i1,1,i2,i3,i4,i5]= 0 #x.f[i1,i2,i3,i4,i5,i6]*ifelse(i6 <= 2,-1,1)
                                y.f[i1,2,i2,i3,i4,i5]= 0
                                y.f[i1,3,i2,i3,i4,i5]= x.f[i1,3,i2,i3,i4,i5]
                                y.f[i1,4,i2,i3,i4,i5]= x.f[i1,4,i2,i3,i4,i5]
                            end
                        end
                    end
                end
            end
        #end

    end

    function mul_1plusγ5x_add!(y::WilsonFermion_4D_nowing_mpi{NC},x::WilsonFermion_4D_nowing_mpi{NC},factor) where {NC}#x = x +(1+gamma_5)/2
        n1,n6,n2,n3,n4,n5 = size(x.f)
        #println("axpby")

        #@inbounds for i6=1:n6
            @inbounds for i5=1:n5
                #it = i5+NDW
                for i4=1:n4
                    #iz = i4+NDW
                    for i3=1:n3
                        #iy = i3+NDW
                        for i2=1:n2
                            #ix = i2+NDW
                            @simd for i1=1:NC

                                #y.f[i1,1,i2,i3,i4,i5]= 0 #x.f[i1,i2,i3,i4,i5,i6]*ifelse(i6 <= 2,-1,1)
                                #y.f[i1,2,i2,i3,i4,i5]= 0
                                y.f[i1,3,i2,i3,i4,i5] += factor*x.f[i1,3,i2,i3,i4,i5]
                                y.f[i1,4,i2,i3,i4,i5] += factor*x.f[i1,4,i2,i3,i4,i5]
                            end
                        end
                    end
                end
            end
        #end

    end

    function mul_1minusγ5x!(y::WilsonFermion_4D_nowing_mpi{NC},x::WilsonFermion_4D_nowing_mpi{NC}) where {NC}#(1-gamma_5)/2
        n1,n6,n2,n3,n4,n5 = size(x.f)
        #println("axpby")

        #@inbounds for i6=1:n6
            @inbounds for i5=1:n5
                #it = i5+NDW
                for i4=1:n4
                    #iz = i4+NDW
                    for i3=1:n3
                        #iy = i3+NDW
                        for i2=1:n2
                            #ix = i2+NDW
                            @simd for i1=1:NC

                                y.f[i1,1,i2,i3,i4,i5]= x.f[i1,1,i2,i3,i4,i5] #x.f[i1,i2,i3,i4,i5,i6]*ifelse(i6 <= 2,-1,1)
                                y.f[i1,2,i2,i3,i4,i5]= x.f[i1,2,i2,i3,i4,i5]
                                y.f[i1,3,i2,i3,i4,i5]= 0#x.f[i1,3,i2,i3,i4,i5]
                                y.f[i1,4,i2,i3,i4,i5]= 0#x.f[i1,3,i2,i3,i4,i5]
                            end
                        end
                    end
                end
            end
        #end

    end

    function mul_1minusγ5x_add!(y::WilsonFermion_4D_nowing_mpi{NC},x::WilsonFermion_4D_nowing_mpi{NC},factor) where {NC}#+(1-gamma_5)/2
        n1,n6,n2,n3,n4,n5 = size(x.f)
        #println("axpby")

        #@inbounds for i6=1:n6
            @inbounds for i5=1:n5
                #it = i5+NDW
                for i4=1:n4
                    #iz = i4+NDW
                    for i3=1:n3
                        #iy = i3+NDW
                        for i2=1:n2
                            #ix = i2+NDW
                            @simd for i1=1:NC

                                #y.f[i1,1,i2,i3,i4,i5]= 0 #x.f[i1,i2,i3,i4,i5,i6]*ifelse(i6 <= 2,-1,1)
                                #y.f[i1,2,i2,i3,i4,i5]= 0
                                y.f[i1,1,i2,i3,i4,i5] += factor*x.f[i1,1,i2,i3,i4,i5]
                                y.f[i1,2,i2,i3,i4,i5] += factor*x.f[i1,2,i2,i3,i4,i5]
                            end
                        end
                    end
                end
            end
        #end
    end

    """
               (       -i )              (       -1 )
     GAMMA1 =  (     -i   )     GAMMA2 = (     +1   )
               (   +i     )              (   +1     )
               ( +i       )              ( -1       )
    """

    function mul_1minusγ1x!(y::WilsonFermion_4D_nowing_mpi{NC},x) where {NC}#(1-gamma_1)
        NX = y.NX
        NY = y.NY
        NZ = y.NZ
        NT = y.NT
        #NC = x.NC
        @inbounds for ic=1:NC
            for it=1:y.PN[4]
                for iz=1:y.PN[3]
                    for iy=1:y.PN[2]
                        @simd for ix=1:y.PN[1]
                            v1 = getvalue(x,ic,1,ix,iy,iz,it) + im*getvalue(x,ic,4,ix,iy,iz,it)
                            v2 = getvalue(x,ic,2,ix,iy,iz,it) + im*getvalue(x,ic,3,ix,iy,iz,it)
                            v3 = getvalue(x,ic,3,ix,iy,iz,it) - im*getvalue(x,ic,2,ix,iy,iz,it)
                            v4 = getvalue(x,ic,4,ix,iy,iz,it) - im*getvalue(x,ic,1,ix,iy,iz,it)
                            setvalue!(y,v1,ic,1,ix,iy,iz,it)#] = v1
                            setvalue!(y,v2,ic,2,ix,iy,iz,it)#] = v2
                            setvalue!(y,v3,ic,3,ix,iy,iz,it)# = v3
                            setvalue!(y,v4,ic,4,ix,iy,iz,it) # = v4
                        end
                    end
                end
            end
        end
    end

    function mul_1plusγ1x!(y::WilsonFermion_4D_nowing_mpi{NC},x) where {NC}#(1+gamma_1)
        NX = y.NX
        NY = y.NY
        NZ = y.NZ
        NT = y.NT
        #NC = x.NC
        @inbounds for ic=1:NC
            for it=1:y.PN[4]
                for iz=1:y.PN[3]
                    for iy=1:y.PN[2]
                        @simd for ix=1:y.PN[1]
                            v1 = getvalue(x,ic,1,ix,iy,iz,it) - im*getvalue(x,ic,4,ix,iy,iz,it)
                            v2 = getvalue(x,ic,2,ix,iy,iz,it) - im*getvalue(x,ic,3,ix,iy,iz,it)
                            v3 = getvalue(x,ic,3,ix,iy,iz,it) + im*getvalue(x,ic,2,ix,iy,iz,it)
                            v4 = getvalue(x,ic,4,ix,iy,iz,it) + im*getvalue(x,ic,1,ix,iy,iz,it)
                            setvalue!(y,v1,ic,1,ix,iy,iz,it)#] = v1
                            setvalue!(y,v2,ic,2,ix,iy,iz,it)#] = v2
                            setvalue!(y,v3,ic,3,ix,iy,iz,it)# = v3
                            setvalue!(y,v4,ic,4,ix,iy,iz,it) # = v4
                        end
                    end
                end
            end
        end
    end

    """
               (       -i )              (       -1 )
     GAMMA1 =  (     -i   )     GAMMA2 = (     +1   )
               (   +i     )              (   +1     )
               ( +i       )              ( -1       )
    """
    function mul_1minusγ2x!(y::WilsonFermion_4D_nowing_mpi{NC},x) where {NC}#(1-gamma_2)
        NX = y.NX
        NY = y.NY
        NZ = y.NZ
        NT = y.NT
        #NC = x.NC
        @inbounds for ic=1:NC
            for it=1:y.PN[4]
                for iz=1:y.PN[3]
                    for iy=1:y.PN[2]
                        @simd for ix=1:y.PN[1]
                            v1 = getvalue(x,ic,1,ix,iy,iz,it) + getvalue(x,ic,4,ix,iy,iz,it)
                            v2 = getvalue(x,ic,2,ix,iy,iz,it) - getvalue(x,ic,3,ix,iy,iz,it)
                            v3 = getvalue(x,ic,3,ix,iy,iz,it) - getvalue(x,ic,2,ix,iy,iz,it)
                            v4 = getvalue(x,ic,4,ix,iy,iz,it) + getvalue(x,ic,1,ix,iy,iz,it)
                            setvalue!(y,v1,ic,1,ix,iy,iz,it)#] = v1
                            setvalue!(y,v2,ic,2,ix,iy,iz,it)#] = v2
                            setvalue!(y,v3,ic,3,ix,iy,iz,it)# = v3
                            setvalue!(y,v4,ic,4,ix,iy,iz,it) # = v4
                        end
                    end
                end
            end
        end
    end

    function mul_1plusγ2x!(y::WilsonFermion_4D_nowing_mpi{NC},x) where {NC}#(1-gamma_5)/2
        NX = y.NX
        NY = y.NY
        NZ = y.NZ
        NT = y.NT
        #NC = x.NC
        @inbounds for ic=1:NC
            for it=1:y.PN[4]
                for iz=1:y.PN[3]
                    for iy=1:y.PN[2]
                        @simd for ix=1:y.PN[1]
                            v1 = getvalue(x,ic,1,ix,iy,iz,it) - getvalue(x,ic,4,ix,iy,iz,it)
                            v2 = getvalue(x,ic,2,ix,iy,iz,it) + getvalue(x,ic,3,ix,iy,iz,it)
                            v3 = getvalue(x,ic,3,ix,iy,iz,it) + getvalue(x,ic,2,ix,iy,iz,it)
                            v4 = getvalue(x,ic,4,ix,iy,iz,it) - getvalue(x,ic,1,ix,iy,iz,it)
                            setvalue!(y,v1,ic,1,ix,iy,iz,it)#] = v1
                            setvalue!(y,v2,ic,2,ix,iy,iz,it)#] = v2
                            setvalue!(y,v3,ic,3,ix,iy,iz,it)# = v3
                            setvalue!(y,v4,ic,4,ix,iy,iz,it) # = v4
                        end
                    end
                end
            end
        end
    end

    """
                   (     -i   )              (     -1   )
         GAMMA3 =  (       +i )     GAMMA4 = (       -1 )
                   ( +i       )              ( -1       )
                   (   -i     )              (   -1     )

    """
    function mul_1minusγ3x!(y::WilsonFermion_4D_nowing_mpi{NC},x) where {NC}#(1-gamma_3)
        NX = y.NX
        NY = y.NY
        NZ = y.NZ
        NT = y.NT
        #NC = x.NC
        @inbounds for ic=1:NC
            for it=1:y.PN[4]
                for iz=1:y.PN[3]
                    for iy=1:y.PN[2]
                        @simd for ix=1:y.PN[1]
                            v1 = getvalue(x,ic,1,ix,iy,iz,it) + im*getvalue(x,ic,3,ix,iy,iz,it)
                            v2 = getvalue(x,ic,2,ix,iy,iz,it) - im*getvalue(x,ic,4,ix,iy,iz,it)
                            v3 = getvalue(x,ic,3,ix,iy,iz,it) - im*getvalue(x,ic,1,ix,iy,iz,it)
                            v4 = getvalue(x,ic,4,ix,iy,iz,it) + im*getvalue(x,ic,2,ix,iy,iz,it)
                            setvalue!(y,v1,ic,1,ix,iy,iz,it)#] = v1
                            setvalue!(y,v2,ic,2,ix,iy,iz,it)#] = v2
                            setvalue!(y,v3,ic,3,ix,iy,iz,it)# = v3
                            setvalue!(y,v4,ic,4,ix,iy,iz,it) # = v4
                        end
                    end
                end
            end
        end
    end

    function mul_1plusγ3x!(y::WilsonFermion_4D_nowing_mpi{NC},x) where {NC}#(1+gamma_3)
        NX = y.NX
        NY = y.NY
        NZ = y.NZ
        NT = y.NT
        #NC = x.NC
        @inbounds for ic=1:NC
            for it=1:y.PN[4]
                for iz=1:y.PN[3]
                    for iy=1:y.PN[2]
                        @simd for ix=1:y.PN[1]
                            v1 = getvalue(x,ic,1,ix,iy,iz,it) - im*getvalue(x,ic,3,ix,iy,iz,it)
                            v2 = getvalue(x,ic,2,ix,iy,iz,it) + im*getvalue(x,ic,4,ix,iy,iz,it)
                            v3 = getvalue(x,ic,3,ix,iy,iz,it) + im*getvalue(x,ic,1,ix,iy,iz,it)
                            v4 = getvalue(x,ic,4,ix,iy,iz,it) - im*getvalue(x,ic,2,ix,iy,iz,it)
                            setvalue!(y,v1,ic,1,ix,iy,iz,it)#] = v1
                            setvalue!(y,v2,ic,2,ix,iy,iz,it)#] = v2
                            setvalue!(y,v3,ic,3,ix,iy,iz,it)# = v3
                            setvalue!(y,v4,ic,4,ix,iy,iz,it) # = v4
                        end
                    end
                end
            end
        end
    end

    """
                   (     -i   )              (     -1   )
         GAMMA3 =  (       +i )     GAMMA4 = (       -1 )
                   ( +i       )              ( -1       )
                   (   -i     )              (   -1     )

    """
    function mul_1minusγ4x!(y::WilsonFermion_4D_nowing_mpi{NC},x) where {NC}#(1-gamma_4)
        NX = y.NX
        NY = y.NY
        NZ = y.NZ
        NT = y.NT
        #NC = x.NC
        @inbounds for ic=1:NC
            for it=1:y.PN[4]
                for iz=1:y.PN[3]
                    for iy=1:y.PN[2]
                        @simd for ix=1:y.PN[1]
                            v1 = getvalue(x,ic,1,ix,iy,iz,it) + getvalue(x,ic,3,ix,iy,iz,it)
                            v2 = getvalue(x,ic,2,ix,iy,iz,it) + getvalue(x,ic,4,ix,iy,iz,it)
                            v3 = getvalue(x,ic,3,ix,iy,iz,it) + getvalue(x,ic,1,ix,iy,iz,it)
                            v4 = getvalue(x,ic,4,ix,iy,iz,it) + getvalue(x,ic,2,ix,iy,iz,it)
                            setvalue!(y,v1,ic,1,ix,iy,iz,it)#] = v1
                            setvalue!(y,v2,ic,2,ix,iy,iz,it)#] = v2
                            setvalue!(y,v3,ic,3,ix,iy,iz,it)# = v3
                            setvalue!(y,v4,ic,4,ix,iy,iz,it) # = v4
                        end
                    end
                end
            end
        end
    end

    function mul_1plusγ4x!(y::WilsonFermion_4D_nowing_mpi{NC},x) where {NC}#(1+gamma_4)
        NX = y.NX
        NY = y.NY
        NZ = y.NZ
        NT = y.NT
        #NC = x.NC
        @inbounds for ic=1:NC
            for it=1:y.PN[4]
                for iz=1:y.PN[3]
                    for iy=1:y.PN[2]
                        @simd for ix=1:y.PN[1]
                            v1 = getvalue(x,ic,1,ix,iy,iz,it) - getvalue(x,ic,3,ix,iy,iz,it)
                            v2 = getvalue(x,ic,2,ix,iy,iz,it) - getvalue(x,ic,4,ix,iy,iz,it)
                            v3 = getvalue(x,ic,3,ix,iy,iz,it) - getvalue(x,ic,1,ix,iy,iz,it)
                            v4 = getvalue(x,ic,4,ix,iy,iz,it) - getvalue(x,ic,2,ix,iy,iz,it)
                            setvalue!(y,v1,ic,1,ix,iy,iz,it)#] = v1
                            setvalue!(y,v2,ic,2,ix,iy,iz,it)#] = v2
                            setvalue!(y,v3,ic,3,ix,iy,iz,it)# = v3
                            setvalue!(y,v4,ic,4,ix,iy,iz,it) # = v4
                        end
                    end
                end
            end
        end
    end


text:
@testset "1658.minimum-operations-to-reduce-x-to-zero.jl" begin
    @test min_operations([1, 1, 14, 2, 3], 5) == 2
    @test min_operations([5, 6, 7, 8, 9], 4) == -1
    @test min_operations([3, 2, 20, 1, 1, 3], 10) == 5
    @test min_operations([8828, 9581, 49, 9818, 9974, 9869, 9991, 10000, 10000, 10000, 9999,
                          9993, 9904, 8819, 1231, 6309], 134365) == 16
end

ext: jl


text:
@testset "Solve getter functions" begin
    ivp = harmonic_oscillator_free()

    alg = Bathe(Δt = 0.1)
    sol = solve(ivp, alg, NSTEPS=100)

    @test dim(sol) == 1

    U = [ui[1] for ui in sol.U]
    @test U == displacements(sol, 1)
    @test_throws ArgumentError displacements(sol, 2)

    V = [vi[1] for vi in sol.U′]
    @test V == velocities(sol, 1)
    @test_throws ArgumentError velocities(sol, 2)

    A = [ai[1] for ai in sol.U′′]
    @test A == accelerations(sol, 1)
    @test_throws ArgumentError accelerations(sol, 2)
end

@testset "Solver options" begin
    ivp = harmonic_oscillator_free()

    alg = Bathe(Δt = 0.1)

    sol = solve(ivp, alg, NSTEPS=100)
    @test times(sol)[end] == 10.0

    sol = solve(ivp, alg, T=10.0)
    @test times(sol)[end] == 10.0

    sol = solve(ivp, alg, tspan=(0.0, 10.0))
    @test times(sol)[end] == 10.0

    sol = solve(ivp, alg, finalTime=10.0)
    @test times(sol)[end] == 10.0

    # test unknown argument name
    @test_throws ArgumentError solve(ivp, alg, final_time=10.0)

    # starting with shifted interval is not implemented
    @test_throws AssertionError solve(ivp, alg, tspan=(1.0, 2.0))
end


text:
module test_2

using AdventOfCode2018.Day2
using Test

@testset "day 2, part 1" begin
	@test parse_input("ab\nbc\ncd\nbr") == ["ab", "bc", "cd", "br"]

	test_input = [
		"abcdef",
		"bababc",
		"abbcde",
		"abcccd",
		"aabcdd",
		"abcdee",
		"ababab",
    ]

	@test Day2.check_id(test_input[1]) == (false, false)
	@test Day2.check_id(test_input[2]) == (true, true)
	@test Day2.check_id(test_input[3]) == (true, false)

	@test part_1(test_input) == 12
end

@testset "day 2, part 2" begin
	test_input = [
		"abcde",
		"fghij",
		"klmno",
		"pqrst",
		"fguij",
		"axcye",
		"wvxyz",
    ]

	@test Day2.is_similar(test_input[1], test_input[2]) == false
	@test Day2.is_similar(test_input[2], test_input[5]) == true
	@test Day2.get_common(test_input[2], test_input[5]) == "fgij"

	@test part_2(test_input) == "fgij"
end
end


text:
# ---
# title: 773. Sliding Puzzle
# id: problem773
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Breadth-first Search
# link: <https://leetcode.com/problems/sliding-puzzle/description/>
# hidden: true
# ---
#
# On a 2x3 `board`, there are 5 tiles represented by the integers 1 through 5,
# and an empty square represented by 0.
#
# A move consists of choosing `0` and a 4-directionally adjacent number and
# swapping it.
#
# The state of the board is _solved_ if and only if the `board` is
# `[[1,2,3],[4,5,0]].`
#
# Given a puzzle board, return the least number of moves required so that the
# state of the board is solved. If it is impossible for the state of the board
# to be solved, return -1.
#
# **Examples:**
#
#
#
#     Input: board = [[1,2,3],[4,0,5]]
#     Output: 1
#     Explanation: Swap the 0 and the 5 in one move.
#
#
#
#     Input: board = [[1,2,3],[5,4,0]]
#     Output: -1
#     Explanation: No number of moves will make the board solved.
#
#
#
#     Input: board = [[4,1,2],[5,0,3]]
#     Output: 5
#     Explanation: 5 is the smallest number of moves that solves the board.
#     An example path:
#     After move 0: [[4,1,2],[5,0,3]]
#     After move 1: [[4,1,2],[0,5,3]]
#     After move 2: [[0,1,2],[4,5,3]]
#     After move 3: [[1,0,2],[4,5,3]]
#     After move 4: [[1,2,0],[4,5,3]]
#     After move 5: [[1,2,3],[4,5,0]]
#
#
#
#     Input: board = [[3,2,4],[1,5,0]]
#     Output: 14
#
#
# **Note:**
#
#   * `board` will be a 2 x 3 array as described above.
#   * `board[i][j]` will be a permutation of `[0, 1, 2, 3, 4, 5]`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module DiskBackedDicts

export DiskBackedDict

using ArgCheck
using JLD2

################################################################################
##### JLD2BlobDict
################################################################################
struct JLD2BlobDict{K,V} <: AbstractDict{K,V}
    path::String
    function JLD2BlobDict{K,V}(path::AbstractString) where {K,V}
        if !(splitext(path)[2] == ".jld2")
            msg = """Path must end with .jld2, got:
            path = $path"""
            throw(ArgumentError(msg))
        end
        new(String(path))
    end
end

function get_dict(obj::JLD2BlobDict{K,V}) where {K,V}
    if !(ispath(obj.path))
        dir = splitdir(obj.path)[1]
        mkpath(dir)
        data = Dict{K,V}()
        JLD2.@save obj.path data
    end
    JLD2.@load(obj.path, data)

    return convert(Dict{K,V},data)
end

function set_dict(obj::JLD2BlobDict{K,V}, data::AbstractDict{K,V}) where {K,V}
    return JLD2.@save obj.path data
end

const PURE_DICT_INTERFACE = [:getindex, :keys, :values, :length, :get, :iterate]
const MUT_DICT_INTERFACE = [:delete!, :setindex!, :get!, :empty!]

for f ∈ PURE_DICT_INTERFACE
    @eval function Base.$f(obj::JLD2BlobDict, args...)
        d = get_dict(obj)
        $f(d, args...)
    end
end

for f in MUT_DICT_INTERFACE
    @eval function Base.$f(obj::JLD2BlobDict, args...)
        d = get_dict(obj)
        ret = $f(d, args...)
        set_dict(obj, d)
        ret
    end
end

function Base.get(f::Base.Callable, o::JLD2BlobDict, key)
    d = get_dict(o)
    get(f,d,key)
end
function Base.get!(f::Base.Callable, o::JLD2BlobDict, key)
    if haskey(o, key)
        nothing
    else
        o[key] = f()
    end
    return o[key]
end

################################################################################
##### FullyCachedDict
################################################################################
"""
    FullyCachedDict([cache,] storage) <: AbstractDict

A fully cached variant of storage. That is `cache` is a copy
of `storage`, but hopefully faster to read. For instance `cache::Dict`
and `storage::JLD2BlobDict`.
"""
struct FullyCachedDict{K,V,C,S} <: AbstractDict{K,V}
    cache::C
    storage::S
    function FullyCachedDict{K,V}(cache::C, storage::S) where {K,V,C,S}
        @argcheck keytype(storage) === keytype(cache) === K
        @argcheck valtype(storage) === valtype(cache) === V
        empty!(cache)
        _merge!(cache, storage)
        return new{K,V,C,S}(cache, storage)
    end
end

function FullyCachedDict(cache, storage)
    @argcheck keytype(storage) === keytype(cache)
    @argcheck valtype(storage) === valtype(cache)
    K = keytype(storage)
    V = valtype(storage)
    return FullyCachedDict{K,V}(cache, storage)
end

function FullyCachedDict(storage::AbstractDict{K,V}) where {K,V}
    cache = Dict{K,V}()
    _merge!(cache, storage)
    return FullyCachedDict{K,V}(cache, storage)
end

for f ∈ PURE_DICT_INTERFACE
    @eval function Base.$f(obj::FullyCachedDict, args...)
        d = obj.cache
        $f(d, args...)
    end
end

for f in MUT_DICT_INTERFACE
    @eval function Base.$f(obj::FullyCachedDict, args...)
        $f(obj.storage, args...)
        $f(obj.cache,   args...)
    end
end

function Base.get(f::Base.Callable, o::FullyCachedDict, key)
    get(f, o.cache, key)
end
function Base.get!(f::Base.Callable, o::FullyCachedDict, key)
    if haskey(o, key)
        val = o[key]
    else
        val = f()
        o[key] = val
    end
    return val
end

################################################################################
##### DiskBackedDict
################################################################################
struct DiskBackedDict{K,V} <: AbstractDict{K,V}
    inner::FullyCachedDict{K,V,Dict{K,V}, JLD2BlobDict{K,V}}
    function DiskBackedDict{K,V}(path::AbstractString) where {K,V}
        storage = JLD2BlobDict{K,V}(path)
        inner = FullyCachedDict(storage)
        return new(inner)
    end
end

DiskBackedDict(path::AbstractString) = DiskBackedDict{Any,Any}(path)

for f in [PURE_DICT_INTERFACE;MUT_DICT_INTERFACE]
    @eval (Base.$f)(d::DiskBackedDict, args...) = $f(d.inner, args...)
end
Base.get(f::Base.Callable, o::DiskBackedDict, key) = Base.get(f, o.inner, key)

################################################################################
##### JLD2FilesStringDict
################################################################################
struct JLD2FilesStringDict{V} <: AbstractDict{String, V}
    root::String
    function JLD2FilesStringDict{V}(path::AbstractString) where {V}
        root = String(path)
        ret = new{V}(root)
        mkpath(valuedir(ret))
        return ret
    end
end

function Base.delete!(o::JLD2FilesStringDict, key)
    rm(valuepath(o, key), force=true)
    return o
end

headerpath(o::JLD2FilesStringDict) = joinpath(o.root, "header.jld2")
valuedir(o::JLD2FilesStringDict) = joinpath(o.root, "values")
valuepath(o::JLD2FilesStringDict, key::AbstractString) = joinpath(valuedir(o), key)

function JLD2FilesStringDict(root::AbstractString)
    V = Any
    JLD2FilesStringDict{V}(root)
end

function Base.empty!(o::JLD2FilesStringDict)
    dir = valuedir(o)
    for filename in readdir(dir)
        path = joinpath(dir, filename)
        rm(path)
    end
    return o
end

function Base.getindex(o::JLD2FilesStringDict, key::AbstractString)
    path = valuepath(o, key)
    ispath(path) || throw(KeyError(key))
    JLD2.@load path value
    return value
end

function Base.keys(o::JLD2FilesStringDict)
    readdir(valuedir(o))
end

function Base.setindex!(o::JLD2FilesStringDict, value, key::AbstractString)
    path = valuepath(o, key)
    value = convert(valtype(o), value)
    JLD2.@save path value
    value
end

Base.length(o::JLD2FilesStringDict) = length(keys(o))

function iterate_pairs_key_based(o)
    ks = keys(o)
    next = iterate(ks)
    if next === nothing
        return nothing
    else
        key, keystate = next
        key => o[key], (keystate=keystate, keys=ks)
    end
end
function iterate_pairs_key_based(o, state)
    next = iterate(state.keys, state.keystate)
    if next === nothing
        return nothing
    else
        key, keystate = next
        key => o[key], (keystate=keystate, keys=state.keys)
    end
end
function getwith_naive!(f, o, key)
    if haskey(o, key)
        return o[key]
    else
        val = convert(valtype(o), f())
        o[key] = val
        return val
    end
end
function getwith_naive(f, o, key)
    if haskey(o, key)
        return o[key]
    else
        f()
    end
end

Base.get!(f::Base.Callable, o::JLD2FilesStringDict, key) = getwith_naive!(f,o,key)
Base.get(f::Base.Callable, o::JLD2FilesStringDict, key) = getwith_naive(f,o,key)

################################################################################
##### JLD2FilesDict
################################################################################
struct JLD2FilesDict{K, V, H} <: AbstractDict{K,V}
    stringdict::JLD2FilesStringDict{Dict{K,V}}
    hash::H
end

Base.get!(f::Base.Callable, o::JLD2FilesDict, key) = getwith_naive!(f,o,key)
Base.get(f::Base.Callable, o::JLD2FilesDict, key) = getwith_naive(f,o,key)

function JLD2FilesDict{K,V}(path::AbstractString, hash=string∘Base.hash) where {K,V}
    H = typeof(hash)
    stringdict = JLD2FilesStringDict{Dict{K,V}}(path)
    return JLD2FilesDict{K,V,H}(stringdict, hash)
end

function JLD2FilesDict(path::AbstractString, hash=string∘Base.hash)
    K = V = Any
    return JLD2FilesDict{K,V}(path, hash)
end

function Base.keys(o::JLD2FilesDict)
    ret = keytype(o)[]
    for skey in keys(o.stringdict)
        append!(ret, keys(o.stringdict[skey]))
    end
    return ret
end

function Base.delete!(o::JLD2FilesDict, key)
    if haskey(o, key)
        skey = _skey(o, key)
        d = o.stringdict[skey]
        if length(d) == 1
            Base.delete!(o.stringdict, skey)
        else
            Base.delete!(d, key)
            o.stringdict[skey] = d
        end
    end
    return o
end
Base.length(o::JLD2FilesDict) = length(keys(o))

Base.iterate(o::Union{JLD2FilesStringDict, JLD2FilesDict}) = iterate_pairs_key_based(o)
Base.iterate(o::Union{JLD2FilesStringDict, JLD2FilesDict}, state) = iterate_pairs_key_based(o, state)

_skey(o, key)::String = o.hash(convert(keytype(o), key))

function Base.get(o::JLD2FilesDict, key, val)
    if haskey(o, key)
        o[key]
    else
        val
    end
end
Base.empty!(o::JLD2FilesDict) = empty!(o.stringdict)

function Base.getindex(o::JLD2FilesDict, key)
    skey = _skey(o, key)
    d = try
        o.stringdict[skey]
    catch err
        if err isa KeyError
            throw(KeyError(key))
        else
            rethrow(err)
        end
    end
    try
        return d[key]
    catch err
        if err isa KeyError
            throw(KeyError(key))
        else
            rethrow(err)
        end
    end
end

function Base.setindex!(o::JLD2FilesDict, val, key)
    K = keytype(o)
    V = valtype(o)
    key = convert(K, key)
    val = convert(V, val)
    skey = _skey(o, key)
    d = get!(o.stringdict, skey) do
        Dict{K,V}()
    end
    d[key] = val
    o.stringdict[skey] = d
    return val
end

################################################################################
##### CachedDict
################################################################################
"""
    CachedDict([cache,] storage) <: AbstractDict

A cached variant of storage. See also [`FullyCachedDict`](@ref).
"""
struct CachedDict{K,V,C,S} <: AbstractDict{K,V}
    cache::C
    storage::S
    function CachedDict{K,V}(cache::C, storage::S) where {K,V,C,S}
        @argcheck keytype(storage) === keytype(cache) === K
        @argcheck valtype(storage) === valtype(cache) === V
        return new{K,V,C,S}(cache, storage)
    end
end

function CachedDict(cache, storage)
    @argcheck keytype(storage) === keytype(cache)
    @argcheck valtype(storage) === valtype(cache)
    K = keytype(storage)
    V = valtype(storage)
    return CachedDict{K,V}(cache, storage)
end

function CachedDict(storage::AbstractDict{K,V}) where {K,V}
    cache = Dict{K,V}()
    return CachedDict{K,V}(cache, storage)
end

function Base.getindex(o::CachedDict, key)
    get!(o.cache, key) do
        o.storage[key]
    end
end
function Base.haskey(o::CachedDict, key)
    haskey(o.cache, key) || haskey(o.storage, key)
end
function Base.setindex!(o::CachedDict, val, key)
    ret = o.storage[key] = val
    o.cache[key] = val
    return ret
end
Base.iterate(o::CachedDict) = iterate_pairs_key_based(o)
Base.iterate(o::CachedDict, state) = iterate_pairs_key_based(o, state)

for f in [:(Base.keys), :(Base.values), :(Base.length)]
    @eval $f(o::CachedDict) = $f(o.storage)
end

function Base.empty!(o::CachedDict)
    empty!(o.cache)
    empty!(o.storage)
end
function Base.delete!(o::CachedDict, key)
    if haskey(o.cache, key)
        Base.delete!(o.cache, key)
    end
    return Base.delete!(o.storage, key)
end
function Base.get(o::CachedDict, key, val)
    get(o.cache, key) do
        get(o.storage, key, val)
    end
end
function Base.get!(o::CachedDict, key, val)
    get!(o.cache, key) do
        get!(o.storage, key, val)
    end
end
function Base.get!(f::Base.Callable, o::CachedDict, key)
    get!(o.cache, key) do
        get!(f, o.storage, key)
    end
end
function Base.get(f::Base.Callable, o::CachedDict, key)
    get(o.cache, key) do
        get(f, o.storage, key)
    end
end

################################################################################
##### merge!
################################################################################

# merge! is an important operation, since it allows to do "batch" otherwise slow
# transactions
function Base.merge!(
        o1::Union{FullyCachedDict,JLD2BlobDict,DiskBackedDict},
        d2::AbstractDict)
    _merge!(o1, d2)
end

function _merge!(d1, d2)
    _merge1!(d1, _batch_reader(d2))
end

_batch_reader(d::JLD2BlobDict) = get_dict(d)
_batch_reader(d::AbstractDict) = d

function _merge1!(o1::JLD2BlobDict, r2)
    d1 = get_dict(o1)
    merge!(d1, r2)
    set_dict(o1, d1)
    o1
end
function _merge1!(d1, r2)
    merge!(d1, r2)
end
function _merge1!(o1::FullyCachedDict, r2)
    _merge1!(o1.cache,   r2)
    _merge1!(o1.storage, r2)
    o1
end
function _merge1!(o1::DiskBackedDict, r2)
    _merge1!(o1.inner, r2)
    o1
end

end # module


text:
using OceanTurb, Printf

@use_pyplot_utils

# # Operations on fields
#
# Generate a grid with 100 points and height H=10. By default the
# base of the grid is at z=-H and the top of the grid is at z=0.
grid = UniformGrid(100, 10)

# Define a function that will generate our field's data.
c0(z) = exp(1.2z)

# Create a field of zeros on the grid.
c = CellField(grid)

# Set c to the function c0.
OceanTurb.set!(c, c0)

# Calculate the derivative of c and set it to cz.
cz = ∂z(c)

# Plot the results
plot(c, label=L"c = e^{1.2z}")
plot(cz, label=L"\partial_z c = 1.2 e^{1.2z}")
legend()

title("A field and it's derivative")
xlabel("\$ c \$ and \$ \\partial_z c \$")
ylabel(L"z")

# # Interpolation
#
# Next we demonstrate how to interpolate from one grid to another,
# respecting the integral budget of a quantity.

grid1 = UniformGrid(30, 10)
grid2 = UniformGrid(10, 10)
grid3 = UniformGrid(3, 10)

c1 = CellField(grid1)
c2 = CellField(grid2)
c3 = CellField(grid3)

OceanTurb.set!(c1, c)
OceanTurb.set!(c2, c)
OceanTurb.set!(c3, c)

@show integral(c)
@show integral(c1)
@show integral(c2)
@show integral(c3)

plot(c, ".", label="high-res \$ c \$")
plot(c1, "^", label="medium-res \$ c \$")
plot(c2, "x", label="low-res \$ c \$")
plot(c3, "+", label="crazy low-res \$ c \$")

legend()

title("\$ e^{1.2z} \$ at different resolutions")
xlabel(L"c")
ylabel(L"z")


text:
function orderedlimits(limits::Union{NamedTuple{(:min, :max), Tuple{Float64, Float64}},
                                     Nothing},
                       limitsname::String)
    if isa(limits, Nothing)
        @info "'$limitsname' limits defined as nothing"
    else
        if limits.max < limits.min
            throw(DataFormatError("$limitsname limits not in ascending order"))
        end
    end

    return limits
end

function getresolution(ts::TimeSeries.TimeArray)

    tstamps = TimeSeries.timestamp(ts)
    timediffs = unique([tstamps[ix]-tstamps[ix-1] for ix in 2:length(tstamps)])

    res = []

    for timediff in timediffs
        if mod(timediff,Dates.Millisecond(Dates.Day(1))) == Dates.Millisecond(0)
            push!(res,Dates.Day(timediff/Dates.Millisecond(Dates.Day(1))))
        elseif mod(timediff,Dates.Millisecond(Dates.Hour(1))) == Dates.Millisecond(0)
            push!(res,Dates.Hour(timediff/Dates.Millisecond(Dates.Hour(1))))
        elseif mod(timediff,Dates.Millisecond(Dates.Minute(1))) == Dates.Millisecond(0)
            push!(res,Dates.Minute(timediff/Dates.Millisecond(Dates.Minute(1))))
        elseif mod(timediff,Dates.Millisecond(Dates.Second(1))) == Dates.Millisecond(0)
            push!(res,Dates.Second(timediff/Dates.Millisecond(Dates.Second(1))))
        else
            throw(DataFormatError("cannot understand the resolution of the timeseries"))
        end
    end

    if length(res) > 1
        throw(DataFormatError(
            "timeseries has non-uniform resolution: this is currently not supported")
        )
    end

    return res[1]
end

"""Throws DataFormatError if the array is not in ascending order."""
function check_ascending_order(array::Array{Int}, name::AbstractString)
    for (i, num) in enumerate(array)
        if i > 1 && num < array[i - 1]
            throw(DataFormatError("$name Numbers are not in ascending order."))
        end
    end

    return
end

"""Checks if a PowerSystemDevice has a field or subfield name."""
function isafield(component::Component, field::Symbol)

    function _wrap(t,d=[])
        fn = fieldnames(typeof(t))
        for n in fn
            push!(d,n)
            f = getfield(t,n)
            if length(fieldnames(typeof(f))) > 0
                _wrap(f,d)
            end
        end
        return d
    end

    allfields = _wrap(component)
    return field in allfields
end


text:
module TestUtil
using LinearAlgebra
using GaussianDistributions: value_gradient_hessian
using Test

@testset "util" begin
    d = 3
    A = randn(d, d)
    A = A'A
    f(x) = dot(x, A, x) / 2
    x = randn(d)
    fx = value_gradient_hessian(f, x)

    @test length(fx) == d^2 + d + 1
    @test fx[1] ≈ f(x)
    @test fx[2:d+1] ≈ A*x
    @test fx[d+2:end] ≈ vec(A)
end
end


text:
abstract type AbstractSignedDistanceFunction end

(::AbstractSignedDistanceFunction)(::Point{3, T}) where {T<:Real} = error("Signed Distance Function for $(typeof(sdf)) has not been defined.")
normal(::AbstractSignedDistanceFunction, ::Point{3, T}) where {T<:Real} = error("Normal of Signed Distance Function for $(typeof(sdf)) has not been defined.")
normalzygote(sdf::AbstractSignedDistanceFunction, point::Point{3, T}) where {T<:Real} =
    Vec3{T}(normalize(Zygote.gradient(sdf, point)[1]))
normalfinitediff(sdf::AbstractSignedDistanceFunction, point::Point{3, T}) where {T<:Real} =
    Vec3{T}(normalize(FiniteDifferences.grad(central_fdm(5, 1), sdf, point)[1]))
normalforwarddiff(sdf::AbstractSignedDistanceFunction, point::Point{3, T}) where {T<:Real} =
    Vec3{T}(normalize(ForwardDiff.gradient(sdf, point)))

struct BoundingVolumeSignedDistanceFunction{T<:Real, BoundingSDF<:AbstractSignedDistanceFunction,
                                         InteriorSDF<:AbstractSignedDistanceFunction} <:
                                         AbstractSignedDistanceFunction
    overshootfactor::T
    boundingSDF::BoundingSDF
    interiorSDF::InteriorSDF
end

overshootfactor(bvsdf::BoundingVolumeSignedDistanceFunction) = bvsdf.overshootfactor
boundingSDF(bvsdf::BoundingVolumeSignedDistanceFunction) = bvsdf.boundingSDF
interiorSDF(bvsdf::BoundingVolumeSignedDistanceFunction) = bvsdf.interiorSDF

function (bvsdf::BoundingVolumeSignedDistanceFunction)(point::Point{3, T}) where {T<:Real}
    bounding_signed_distance = boundingSDF(bvsdf)(point)
    if bounding_signed_distance > T(overshootfactor(bvsdf))
        bounding_signed_distance
    else
        interiorSDF(bvsdf)(point)
    end
end

function normal(bvsdf::BoundingVolumeSignedDistanceFunction, point::Point{3, T}) where {T<:Real}
    bounding_signed_distance = boundingSDF(bvsdf)(point)
    if bounding_signed_distance > T(overshootfactor(bvsdf))
        normal(boundingSDF(bvsdf), point)
    else
        normal(interiorSDF(bvsdf), point)
    end
end

struct UnionSignedDistanceFunction{T<:Tuple{Vararg{AbstractSignedDistanceFunction}}} <: AbstractSignedDistanceFunction
    SDFs::T
end

SDFs(union_sdf::UnionSignedDistanceFunction) = union_sdf.SDFs

function (union_sdf::UnionSignedDistanceFunction)(point::Point{3, T}) where {T<:Real}
    minimum(sdf(point) for sdf in SDFs(union_sdf))
end

function normal(union_sdf::UnionSignedDistanceFunction, point::Point{3, T}) where {T<:Real}
    min_sdf_index = argmin(sdf(point) for sdf in SDFs(union_sdf))
    normal(SDFs(union_sdf)[min_sdf_index], point)
end

struct IntersectSignedDistanceFunction{T<:Tuple{Vararg{AbstractSignedDistanceFunction}}} <: AbstractSignedDistanceFunction
    SDFs::T
end

SDFs(intersect_sdf::IntersectSignedDistanceFunction) = intersect_sdf.SDFs

function (intersect_sdf::IntersectSignedDistanceFunction)(point::Point{3, T}) where {T<:Real}
    maximum(sdf(point) for sdf in SDFs(intersect_sdf))
end

function normal(intersect_sdf::IntersectSignedDistanceFunction, point::Point{3, T}) where {T<:Real}
    max_sdf_index = argmax(sdf(point) for sdf in SDFs(intersect_sdf))
    normal(SDFs(intersect_sdf)[max_sdf_index], point)
end

struct NegatedSignedDistanceFunction{T<:AbstractSignedDistanceFunction} <: AbstractSignedDistanceFunction
    sdf::T
end

sdf(negated_sdf::NegatedSignedDistanceFunction) = negated_sdf.sdf

function (negated_sdf::NegatedSignedDistanceFunction)(point::Point{3, T}) where {T<:Real}
    -(sdf(negated_sdf)(point))
end

function normal(negated_sdf::NegatedSignedDistanceFunction, point::Point{3, T}) where {T<:Real}
    -(normal(sdf(negated_sdf), point))
end

struct TransformedSignedDistanceFunction{Transform<:Transformation, SDF<:AbstractSignedDistanceFunction} <: AbstractSignedDistanceFunction
    transform::Transform
    invtransform::Transform
    sdf::SDF
end

transform(transformed_sdf::TransformedSignedDistanceFunction) = transformed_sdf.transform
invtransform(transformed_sdf::TransformedSignedDistanceFunction) = transformed_sdf.invtransform
sdf(transformed_sdf::TransformedSignedDistanceFunction) = transformed_sdf.sdf
transform_sdf(transformation::Transformation, sdf::AbstractSignedDistanceFunction) = TransformedSignedDistanceFunction(transformation, inv(transformation), sdf)
function transform_sdf(transformation::Transformation, transformed_sdf::TransformedSignedDistanceFunction)
    # collapse composed transformations
    composed_transformation = compose(transformation, transform(transformed_sdf))
    TransformedSignedDistanceFunction(composed_transformation, inv(composed_transformation), sdf(transformed_sdf))
end

(transformed_sdf::TransformedSignedDistanceFunction)(point::Point{3, T}) where {T<:Real} = sdf(transformed_sdf)(invtransform(transformed_sdf)(point))

function normal(transformed_sdf::TransformedSignedDistanceFunction, point::Point{3, T}) where {T<:Real}
    Vec{3, T}(normalize(LinearMap(transform_deriv(invtransform(transformed_sdf), invtransform(transformed_sdf)(point)))(invtransform(transformed_sdf)(point))))
end


text:
"""
$(TYPEDEF)
"""
struct StandardDDEProblem end

"""
$(TYPEDEF)
"""
struct DDEProblem{uType,tType,lType,lType2,isinplace,P,F,H,K,PT} <:
                          AbstractDDEProblem{uType,tType,lType,isinplace}
  f::F
  u0::uType
  h::H
  tspan::tType
  p::P
  constant_lags::lType
  dependent_lags::lType2
  kwargs::K
  neutral::Bool
  order_discontinuity_t0::Int
  problem_type::PT

  @add_kwonly function DDEProblem{iip}(f::AbstractDDEFunction{iip}, u0, h, tspan, p=NullParameters();
                                       constant_lags = (),
                                       dependent_lags = (),
                                       neutral = f.mass_matrix !== I && det(f.mass_matrix) != 1,
                                       order_discontinuity_t0 = 0,
                                       problem_type = StandardDDEProblem(),
                                       kwargs...) where {iip}
    _tspan = promote_tspan(tspan)
    new{typeof(u0),typeof(_tspan),typeof(constant_lags),typeof(dependent_lags),isinplace(f),
        typeof(p),typeof(f),typeof(h),typeof(kwargs),typeof(problem_type)}(
          f, u0, h, _tspan, p, constant_lags, dependent_lags, kwargs, neutral,
          order_discontinuity_t0, problem_type)
  end

  function DDEProblem{iip}(f::AbstractDDEFunction{iip}, h, tspan::Tuple, p=NullParameters();
                           order_discontinuity_t0 = 1, kwargs...) where iip
    DDEProblem{iip}(f, h(p, first(tspan)), h, tspan, p;
                    order_discontinuity_t0 = max(1, order_discontinuity_t0), kwargs...)
  end

  function DDEProblem{iip}(f, args...; kwargs...) where iip
    DDEProblem{iip}(convert(DDEFunction{iip}, f), args...; kwargs...)
  end
end

DDEProblem(f, args...; kwargs...) =
  DDEProblem(convert(DDEFunction, f), args...; kwargs...)

DDEProblem(f::AbstractDDEFunction, args...; kwargs...) =
  DDEProblem{isinplace(f)}(f, args...; kwargs...)

"""
$(TYPEDEF)
"""
abstract type AbstractDynamicalDDEProblem end

"""
$(TYPEDEF)
"""
struct DynamicalDDEProblem{iip} <: AbstractDynamicalDDEProblem end

# u' = f1(v,h)
# v' = f2(t,u,h)
"""
    DynamicalDDEProblem(f::DynamicalDDEFunction,v0,u0,tspan,p=NullParameters(),callback=CallbackSet())

Define a dynamical DDE problem from a [`DynamicalDDEFunction`](@ref).
"""
function DynamicalDDEProblem(f::DynamicalDDEFunction,v0,u0,h,tspan,p=NullParameters();dependent_lags=(),kwargs...)
  DDEProblem(f,ArrayPartition(v0,u0),h,tspan,p;
             problem_type=DynamicalDDEProblem{isinplace(f)}(),
             dependent_lags=ntuple(i->(u,p,t)->dependent_lags[i](u[1],u[2],p,t),length(dependent_lags)),
             kwargs...)
end
function DynamicalDDEProblem(f::DynamicalDDEFunction,h,tspan,p=NullParameters();kwargs...)
  DynamicalDDEProblem(f,h(p,first(tspan))...,h,tspan,p;kwargs...)
end
function DynamicalDDEProblem(f1,f2,args...;kwargs...)
  DynamicalDDEProblem(DynamicalDDEFunction(f1,f2),args...;kwargs...)
end

"""
    DynamicalDDEProblem{isinplace}(f1,f2,v0,u0,h,tspan,p=NullParameters(),callback=CallbackSet())

Define a dynamical DDE problem from the two functions `f1` and `f2`.

# Arguments
* `f1` and `f2`: The functions in the DDE.
* `v0` and `u0`: The initial conditions.
* `h`: The initial history function.
* `tspan`: The timespan for the problem.
* `p`: Parameter values for `f1` and `f2`.
* `callback`: A callback to be applied to every solver which uses the problem. Defaults to nothing.

`isinplace` optionally sets whether the function is inplace or not.
This is determined automatically, but not inferred.
"""
function DynamicalDDEProblem{iip}(f1,f2,args...;kwargs...) where iip
  DynamicalDDEProblem(DynamicalDDEFunction{iip}(f1,f2),args...;kwargs...)
end

# u'' = f(du,u,h,p,t)
"""
$(TYPEDEF)
"""
struct SecondOrderDDEProblem{iip} <: AbstractDynamicalDDEProblem end
function SecondOrderDDEProblem(f,args...;kwargs...)
  iip = isinplace(f,6)
  SecondOrderDDEProblem{iip}(f,args...;kwargs...)
end

"""
    SecondOrderDDEProblem{isinplace}(f,du0,u0,h,tspan,p=NullParameters(),callback=CallbackSet())

Define a second order DDE problem with the specified function.

# Arguments
* `f`: The function for the second derivative.
* `du0`: The initial derivative.
* `u0`: The initial condition.
* `h`: The initial history function.
* `tspan`: The timespan for the problem.
* `p`: Parameter values for `f`.
* `callback`: A callback to be applied to every solver which uses the problem. Defaults to nothing.

`isinplace` optionally sets whether the function is inplace or not.
This is determined automatically, but not inferred.
"""
function SecondOrderDDEProblem{iip}(f,args...;kwargs...) where iip
  if iip
    f2 = function (du,v,u,h,p,t)
      du .= v
    end
  else
    f2 = function (v,u,h,p,t)
      v
    end
  end
  DynamicalDDEProblem{iip}(f,f2,args...;problem_type=SecondOrderDDEProblem{iip}(),kwargs...)
end
function SecondOrderDDEProblem(f::DynamicalDDEFunction,args...;kwargs...)
  iip = isinplace(f.f1, 6)
  if f.f2.f === nothing
    if iip
      f2 = function (du,v,u,h,p,t)
        du .= v
      end
    else
      f2 = function (v,u,h,p,t)
        v
      end
    end
    return DynamicalDDEProblem(DynamicalDDEFunction{iip}(f.f1,f2;mass_matrix=f.mass_matrix,analytic=f.analytic),
                               args...;problem_type=SecondOrderDDEProblem{iip}(),kwargs...)
  else
    return DynamicalDDEProblem(DynamicalDDEFunction{iip}(f.f1,f.f2;mass_matrix=f.mass_matrix,analytic=f.analytic),
                               args...;problem_type=SecondOrderDDEProblem{iip}(),kwargs...)
  end
end


text:
using BenchmarkTools
using Aligator

include("singlepath.jl")
include("multipath.jl")

singlepath = [:cohencu, :freire1, :freire2, :(petter(1)), :(petter(2)), :(petter(3)), :(petter(4))]
multipath = [:divbin, :euclidex, :fermat, :knuth, :lcm, :mannadiv, :wensley]

const ijcar18 = BenchmarkGroup()
ijcar18["singlepath"] = BenchmarkGroup()
ijcar18["multipath"] = BenchmarkGroup()

macro createbenchmarks(suite, instances)
    for loop in eval(instances)
        loopstr = eval(loop)
        suite = eval(suite)
        suite[string(loop)] = @benchmarkable aligator($(loopstr))
    end
end

@createbenchmarks ijcar18["singlepath"] singlepath
@createbenchmarks ijcar18["multipath"] multipath

tune!(ijcar18)
results = run(ijcar18)

showall(results)

# results["singlepath"]
# results["multipath"]