text:
using SymPy

using Oceananigans.Advection:
    eno_coefficients, optimal_weno_weights, β,
    weno_flux_x, weno_flux_y, weno_flux_z,
    left_biased_interpolate_xᶠᵃᵃ, left_biased_interpolate_yᵃᶠᵃ, left_biased_interpolate_zᵃᵃᶠ

function print_eno_interpolant(k, r)
    cs = eno_coefficients(k, r)

    ssign(n) = n >= 0 ? "+" : "-"
    ssubscript(n) = n == 0 ? "i"  : "i" * ssign(n) * string(abs(n))

    print("u(i+1//2) = ")
    for (j, c) in enumerate(cs)
        c_s = ssign(c) * " " * string(abs(c))
        ss_s = ssubscript(j-r-1)
        print(c_s * " u(" * ss_s * ") ")
    end
    print("\n")
end

function print_eno_interpolants(k)
    for r in -1:k-1
        print_eno_interpolant(k, r)
    end
    Γ = optimal_weno_weights(k)
    println("Optimal weights γ: $Γ")
end

@testset "WENO reconstruction" begin
    @info "Testing WENO reconstruction..."

    # Janky regression test

    weno5 = WENO(5)

    N = 10
    a = collect(Float64, (1:N) .^ 2)

    ax = reshape(a, (N, 1, 1))
    ay = reshape(a, (1, N, 1))
    az = reshape(a, (1, 1, N))

    correct = [73//6, 121//6, 181//6, 253//6]

    inds = 4:N-3
    @test all([weno_flux_x(i, 1, 1, weno5, ax) for i in inds] .≈ correct)
    @test all([weno_flux_y(1, j, 1, weno5, ay) for j in inds] .≈ correct)
    @test all([weno_flux_z(1, 1, k, weno5, az) for k in inds] .≈ correct)

    @test rationalize.([left_biased_interpolate_xᶠᵃᵃ(i, 1, 1, nothing, WENO5(), ax) for i in inds]) == correct
    @test rationalize.([left_biased_interpolate_yᵃᶠᵃ(1, j, 1, nothing, WENO5(), ay) for j in inds]) == correct
    @test rationalize.([left_biased_interpolate_zᵃᵃᶠ(1, 1, k, nothing, WENO5(), az) for k in inds]) == correct

    @testset "ENO reconstruction weights" begin
        @testset "WENO-3" begin
            @info "WENO-3 coefficients [Compare with Table 2.1 of Shu (1998)]:"
            print_eno_interpolants(2)

            @test eno_coefficients(2, -1) == [ 3//2, -1//2]
            @test eno_coefficients(2,  0) == [ 1//2,  1//2]
            @test eno_coefficients(2,  1) == [-1//2,  3//2]
        end

        @testset "WENO-5" begin
            @info "WENO-5 coefficients [Compare with Table 2.1 of Shu (1998) and equation (2.15) from Shu (2009)]:"
            print_eno_interpolants(3)

            @test eno_coefficients(3, -1) == [11//6, -7//6,  1//3]
            @test eno_coefficients(3,  0) == [ 1//3,  5//6, -1//6]
            @test eno_coefficients(3,  1) == [-1//6,  5//6,  1//3]
            @test eno_coefficients(3,  2) == [ 1//3, -7//6, 11//6]
        end

        @testset "WENO-7" begin
            @info "WENO-7 coefficients [Compare with Table 2.1 of Shu (1998)]:"
            print_eno_interpolants(4)

            @test eno_coefficients(4, -1) == [25//12, -23//12,  13//12, -1//4 ]
            @test eno_coefficients(4,  0) == [ 1//4,   13//12,  -5//12,  1//12]
            @test eno_coefficients(4,  1) == [-1//12,   7//12,   7//12, -1//12]
            @test eno_coefficients(4,  2) == [ 1//12,  -5//12,  13//12,  1//4 ]
            @test eno_coefficients(4,  3) == [-1//4,   13//12, -23//12, 25//12]
        end

        @testset "WENO-9" begin
            @info "WENO-9 coefficients [Compare with Table 2.1 of Shu (1998) and equation (2.14) of Shu (2009)]:"
            print_eno_interpolants(5)

            @test eno_coefficients(5, -1) == [ 137//60, -163//60, 137//60,  -21//20,   1//5 ]
            @test eno_coefficients(5,  0) == [   1//5,    77//60, -43//60,   17//60,  -1//20]
            @test eno_coefficients(5,  1) == [  -1//20,    9//20,  47//60,  -13//60,   1//30]
            @test eno_coefficients(5,  2) == [   1//30,  -13//60,  47//60,    9//20,  -1//20]
            @test eno_coefficients(5,  3) == [  -1//20,   17//60, -43//60,   77//60,   1//5 ]
            @test eno_coefficients(5,  4) == [   1//5,   -21//20, 137//60, -163//60, 137//60]
        end
    end

    @testset "WENO optimal weights" begin
        @testset "WENO-3" begin
            # Compare with Table II of Jiang & Shu (1996).
            @test optimal_weno_weights(2) == [2//3, 1//3]
        end

        @testset "WENO-5" begin
            # Compare with equation (2.15) of Shu (2009).
            @test optimal_weno_weights(3) == [3//10, 3//5, 1//10]
        end

        @testset "WENO-7" begin end

        @testset "WENO-9" begin
            @test optimal_weno_weights(5) == [5//126, 20//63, 10//21, 10//63, 1//126]
        end
    end

    @testset "WENO smoothness indicators" begin
        @testset "WENO-5" begin
            @vars ϕᵢ₋₂  ϕᵢ₋₁ ϕᵢ ϕᵢ₊₁ ϕᵢ₊₂

            β₀_correct = 13//12 * (ϕᵢ₋₂ - 2ϕᵢ₋₁ + ϕᵢ  )^2 + 1//4 * ( ϕᵢ₋₂ - 4ϕᵢ₋₁ + 3ϕᵢ  )^2 |> expand
            β₁_correct = 13//12 * (ϕᵢ₋₁ - 2ϕᵢ   + ϕᵢ₊₁)^2 + 1//4 * ( ϕᵢ₋₁         -  ϕᵢ₊₁)^2 |> expand
            β₂_correct = 13//12 * (ϕᵢ   - 2ϕᵢ₊₁ + ϕᵢ₊₂)^2 + 1//4 * (3ϕᵢ   - 4ϕᵢ₊₁ +  ϕᵢ₊₂)^2 |> expand

            @test β(3, 0, (ϕᵢ, ϕᵢ₋₁, ϕᵢ₋₂)) |> expand == β₀_correct
            @test β(3, 1, (ϕᵢ₊₁, ϕᵢ, ϕᵢ₋₁)) |> expand == β₁_correct
            @test β(3, 2, (ϕᵢ₊₂, ϕᵢ₊₁, ϕᵢ)) |> expand == β₂_correct
        end
    end
end


text:
@testset "traits.jl" begin
    for t in (OneToOne(), OneToMany(), ManyToOne(), ManyToMany())
        @test t isa FeatureTransforms.Cardinality
    end
end

ext: jl


text:
using FastGaussQuadrature

#=
Gaußian quadrature:

\[\int\limits_a^b dx\;f(x)\approx
\frac{b-a}{2}\sum_{i=1}^n w_i
f\left(\frac{b-a}{2}x_i+\frac{a+b}{2}\right)\]
=#
function lgwt!(x::AbstractVector{T}, w::AbstractVector{T},
               xs::AbstractVector{T}, ws::AbstractVector{T},
               i,
               a::T=zero(T), b::T=one(T)) where {T<:Real}
    n = length(x)
    interval = (i-1)*n+1 : i*n
    xs[interval] = 0.5((b-a)*x .+ a .+ b)
    ws[interval] = 0.5(b-a)*w
    nothing
end

function lgwt(t::AbstractKnotSet{T}) where T
    k = order(t)
    N = numintervals(t)
    x, w = gausslegendre(k+1)
    xo = zeros(T, N*length(x))
    wo = zeros(T, N*length(x))

    for i = 1:N
        lgwt!(x, w,
              xo, wo, i,
              t[k+i-1], t[k+i])
    end

    xo,wo
end


text:
# BSE kernel for the t channel
function compute_t_BSE!(
    Λ    :: Float64,
    buff :: Matrix{Float64},
    v    :: Float64,
    dv   :: Float64,
    t    :: Float64,
    vt   :: Float64,
    vtp  :: Float64,
    r    :: Reduced_lattice,
    m    :: Mesh,
    a    :: Action_su2,
    temp :: Array{Float64, 3}
    )    :: Nothing

    # get propagator and overlap
    p       = -get_propagator(Λ, v + 0.5 * t, v - 0.5 * t, m, a)
    overlap = r.overlap

    # get buffers for left non-local vertex
    bs1 = get_buffer_s( v + vt, 0.5 * (-t - v + vt), 0.5 * (-t + v - vt), m)
    bt1 = get_buffer_empty()
    bu1 = get_buffer_u(-v + vt, 0.5 * (-t + v + vt), 0.5 * ( t + v + vt), m)

    # get buffers for right non-local vertex
    bs2 = get_buffer_s(v + vtp, 0.5 * (-t + v - vtp), 0.5 * (-t - v + vtp), m)
    bt2 = get_buffer_t(t, v, vtp, m)
    bu2 = get_buffer_u(v - vtp, 0.5 * (-t + v + vtp), 0.5 * ( t + v + vtp), m)

    # get buffers for local left vertex
    bs3 = get_buffer_s( v + vt, 0.5 * (t + v - vt), 0.5 * (-t + v - vt), m)
    bt3 = get_buffer_t(-v + vt, 0.5 * (t + v + vt), 0.5 * (-t + v + vt), m)
    bu3 = get_buffer_empty()

    # get buffers for local right vertex
    bs4 = get_buffer_s(v + vtp, 0.5 * (t - v + vtp), 0.5 * (-t - v + vtp), m)
    bt4 = get_buffer_t(v - vtp, 0.5 * (t + v + vtp), 0.5 * (-t + v + vtp), m)
    bu4 = get_buffer_u(t, vtp, v, m)

    # cache local vertex values
    v3s, v3d = get_Γ(1, bs3, bt3, bu3, r, a, ch_u = false)
    v4s, v4d = get_Γ(1, bs4, bt4, bu4, r, a)

    # cache vertex values for all lattice sites in temporary buffer
    get_Γ_avx!(r, bs1, bt1, bu1, a, temp, 1, ch_t = false)
    get_Γ_avx!(r, bs2, bt2, bu2, a, temp, 2)

    # compute contributions for all lattice sites
    for i in eachindex(r.sites)
        # read cached values for site i
        v1s = temp[i, 1, 1]; v1d = temp[i, 2, 1]
        v2s = temp[i, 1, 2]; v2d = temp[i, 2, 2]

        # compute contribution at site i
        Γs = -p * (-1.0 * v1s * v4s + v1s * v4d - 1.0 * v3s * v2s + v3d * v2s)
        Γd = -p * ( 3.0 * v1d * v4s + v1d * v4d + 3.0 * v3s * v2d + v3d * v2d)

        # determine overlap for site i
        overlap_i = overlap[i]

        # determine range for inner sum
        Range = size(overlap_i, 1)

        # compute inner sum
        @turbo unroll = 1 for j in 1 : Range
            # read cached values for inner site
            v1s = temp[overlap_i[j, 1], 1, 1]; v1d = temp[overlap_i[j, 1], 2, 1]
            v2s = temp[overlap_i[j, 2], 1, 2]; v2d = temp[overlap_i[j, 2], 2, 2]

            # compute contribution at inner site
            Γs += -p * (-2.0) * overlap_i[j, 3] * (2.0 * a.S) * v1s * v2s
            Γd += -p * (-2.0) * overlap_i[j, 3] * (2.0 * a.S) * v1d * v2d
        end

        # parse result to output buffer
        buff[1, i] += dv * Γs
        buff[2, i] += dv * Γd
    end

    return nothing
end


text:
using Bolt
using PyPlot
using Unitful, UnitfulAstro, NaturallyUnitful
using Interpolations

par = CosmoParams()
bg = Background(par)
ih = IonizationHistory(Peebles(), par, bg)
x_grid = bg.x_grid

clf()
fig, ax = subplots(1,2,figsize=(10,5))
ax[1].plot(x_grid, ih.τ.(x_grid), "-", label=raw"$\tau$")
ax[1].plot(x_grid, abs.(ih.τ′.(x_grid)), "--", label=raw"$|\tau^\prime|$")
ax[2].plot(x_grid, ih.g̃.(x_grid), "-", label=raw"$\tilde{g}$")
ax[2].plot(x_grid, ih.g̃′.(x_grid) ./ 10, "--", label=raw"$\tilde{g}\prime/10$")
ax[2].plot(x_grid, ih.g̃′′.(x_grid) ./ 300, "--", label=raw"$\tilde{g}\prime/300$")
ax[1].set_yscale("log")
ax[1].legend()
ax[1].set_xlabel(raw"$x$")
ax[2].set_xlim(-8.0, -6.0)
ax[2].set_ylim(-3.5, 5.5)
ax[2].legend()
ax[2].set_xlabel(raw"$x$")
tight_layout()
gcf()


text:
module List
export List!, Cons, Nil
using MLStyle
using MLStyle.Data
import Base: ^, map, foreach, iterate, isempty, getindex

@data List!{T} begin
    Nil{T}
    Cons{T}(head :: T, tail :: List!{T})
end

# NO SFINAE in Julia...

Base.convert(::Type{List!{G}}, it :: Nil{T}) where {G, T <: G} = Nil{G}()

Base.convert(::Type{List!{G}}, it :: Cons{T}) where {G, T <: G} =
    Cons{G}(
        convert(G, it.head),
        convert(List!{G}, it.tail))

Cons(head :: T, tail :: List!{G}) where {G, T} =
    Cons{typejoin(T, G)}(head, tail)

function (^)(head :: T, tail :: List!{G}) where {G, T}
    Cons{typejoin(T, G)}(head, tail)
end

function List!(args :: Vararg{T, N}) where {T, N}
    ret = Nil{T}()
    for i = 0:N-1
        ret = Cons{T}(args[end - i], ret)
    end
    ret
end

function (::Type{List!{T1}})(args :: Vararg{T2, N}) where {T1, T2, N}
    T =
        if T1 === Any && T2 !== Any
            T2
        else
            T1
        end

    ret = Nil{T}()

    for i = 0:N-1
        ret = Cons{T}(args[end - i], ret)
    end
    ret
end

List!(args...) =
    begin
        ret = Nil{Any}()
        for i = 0:length(args)-1
            ret = Cons{Any}(args[end - i], ret)
        end
        ret
    end

Nil() = Nil{Any}()

Base.iterate(lst :: List!{T}) where T =
    iterate(lst, lst)

function Base.iterate(::List!{T}, lst :: List!{T})::Union{Nothing, Tuple{T, List!{T}}} where T
    @match lst begin
        Nil{T}()      => nothing
        Cons{T}(a, b) => (a, b)
    end
end

function Base.map(fn :: Fun{T, R}, lst :: List!{T}) :: List!{R} where {T, R}
    ret = Nil{R}()
    for each in lst
       ret = Cons{R}(fn(each), ret)
    end
    ret
end

function Base.map(fn :: Function, lst :: List!{T}) :: List! where T

    ret = Nil()
    for each in lst
        ret = Cons(fn(each), ret)
    end
    ret

end

function Base.foreach(fn :: Fun{T, R}, lst :: List!{T}) :: Nothing where {T, R}
    for each in lst
        fn(each)
    end
end

function Base.foreach(fn :: Function, lst :: List!{T}) :: Nothing  where T
    for each in lst
        fn(each)
    end
end

function _flatten_cons_marker(ast :: Expr, res :: Vector)
    if ast.head === :call && ast.args[1] === :^
        push!(res, ast.args[2])
        _flatten_cons_marker(ast.args[3], res)
    else
        push!(res, ast)
    end
end

function _flatten_cons_marker(ast :: Symbol, res :: Vector)
    push!(res, ast)
end

function _flatten_cons_marker(ast)
    if isa(ast, Expr)
        res = []
        _flatten_cons_marker(ast, res)
    else
        [ast]
    end
end

function _match(args, guard, tag, mod)
    l, r = args

    let l    = pattern_match(l, nothing, :($tag.head), mod),
        tail = :($tag.tail),
        r    = if isa(r, Expr) && r.head == :call && r.args[1] === :^
                   _match(r.args[2:end], nothing, tail, mod)
               else
                   pattern_match(r, nothing, tail, mod)
               end

        :($l && $r)

    end |>
    function (last)
        if guard === nothing
            last
        else
            :($last && $guard)
        end
    end
end

PatternDef.App(^) do args, guard, tag, mod
    check_ty = :($isa($tag, $Cons))
    _match(args, nothing, tag, mod) |>
    function (last)
        if guard === nothing
            last
        else
            :($last && $guard)
        end
    end
end

# TODO:
# Base.reduce, Base.flatmap, Base.all, Base.any, Base.isempty, Base.copy
# head, tail, find, take, zip...

function isempty(lst :: Cons{T}) where T
    false
end

function isempty(lst :: Nil{T}) where T
    true
end

function head(lst :: Cons{T}) :: T where T
    lst.head
end

function head(lst :: Nil{T}) where T
    ArgumentError("Cannot get head from zero sized list.") |> throw
end

function try_head(lst :: Cons{T}) where T
    Some{T}(lst.head)
end

function try_head(lst :: Nil{T}) where T
    nothing
end

function tail(lst :: Cons{T}) where T
    lst.tail
end

function tail(lst :: Nil{T}) where T
    ArgumentError("Cannot get tail from zero sized list.") |> throw
end

function try_tail(lst :: Cons{T}) where T
    Some{T}(lst.tail)
end

function try_tail(lst :: Nil{T}) where T
    nothing
end

end


text:
include("classes.jl")
include("private/private.jl")

# COMPONENTS
struct Transform
    m_position::AbstractArray{Float32}
    m_rotation::AbstractArray{Float32}
    m_scale::AbstractArray{Float32}
    model_matrix::AbstractArray{AbstractArray{Float32}}
    m_right_axis::AbstractArray{Float32}
	m_front_axis::AbstractArray{Float32}
	m_up_axis::AbstractArray{Float32}

    Transform() = new([0, 0, 0], [0, 0, 0], [0, 0, 0],
    [[1, 0, 0, 0],
     [0, 1, 0, 0],
     [0, 0, 1, 0],
     [0, 0, 0, 1]], [1, 0, 0], [0, 0, -1], [0, 1, 0])
end

# FIX THE type of ptr members
struct SpriteRenderer
    m_texture::Texture #ptr
    m_vertices::Vertex #ptr
    m_color::AbstractArray{Float32}
    m_material::Material #ptr
    m_shader_program::_private.J3DShaderProgram
    m_viabo::_private.VertexIndexArrayBufferObject

    # FIX THE type of ptr members
    SpriteRenderer() = new(Texture(), Vertex(), [0, 0, 0, 1],
    Material(), _private.J3DShaderProgram(), _private.VertexIndexArrayBufferObject())
end

struct Camera

end


text:
using GridInterpolations
using LocalFunctionApproximation
using LocalApproximationValueIteration
using POMDPs
using POMDPModels
using POMDPModelTools
using StaticArrays
using JLD2, FileIO
using Random
using Logging
using Distributions
using CMSSPs

# Read parameters from file
params_fn = ARGS[1]
inhor_fn = ARGS[2]
outhor_fn = ARGS[3]

rng = MersenneTwister(10)

# Create parameters and then CMSSP
params = toy2d_parse_params(params_fn)
cmssp = create_toy2d_cmssp(params)

# Set up the modal mdp
# Need to do it for just one mode
modal_mdp = Toy2DModalMDPType(1, params, cmssp.control_actions, 0.8, params.horizon_limit)

# Create grid interpolator
xy_spacing = polyspace_symmetric(1.0, params.axis_vals, 3)
xy_grid = RectangleGrid(xy_spacing, xy_spacing)
lfa = LocalGIFunctionApproximator(xy_grid)

# Compute and update the terminal cost
compute_terminalcost_localapprox!(modal_mdp, cmssp, 1, lfa)

# Compute policy with localApproxVI
modal_policy = finite_horizon_VI_localapprox!(modal_mdp, lfa, true,
                                              params.num_generative_samples, rng)

# Compute min_cost_per_horizon
compute_min_value_per_horizon_localapprox!(modal_policy)

# Save horizon policy to file
save_modal_horizon_policy_localapprox(modal_policy, inhor_fn,outhor_fn, modal_mdp)


text:
using Test

using JacobiDavidson, LinearAlgebra, LinearMaps

@testset "Schur permutations" begin
    A = rand(ComplexF64, 10, 10)
    B = rand(ComplexF64, 10, 10)
    F = schur(A, B)

    τ = 1.0 + 0.0im
    eigenvalues = F.alpha ./ F.beta

    # Test whether the ordering is correct
    closest_to_one = first(schur_permutation(Near(τ), F, 1:1))
    @test abs.(eigenvalues[closest_to_one] - τ) ≤ minimum(abs.(eigenvalues .- τ))

    largest_magnitude = first(schur_permutation(LargestMagnitude(τ), F, 1:1))
    @test abs(eigenvalues[largest_magnitude]) ≥ maximum(abs.(eigenvalues))

    smallest_magnitude = first(schur_permutation(SmallestMagnitude(τ), F, 1:1))
    @test abs(eigenvalues[smallest_magnitude]) ≤ minimum(abs.(eigenvalues))

    largest_real_part = first(schur_permutation(LargestRealPart(τ), F, 1:1))
    @test real(eigenvalues[largest_real_part]) ≥ maximum(real.(eigenvalues))

    smallest_real_part = first(schur_permutation(SmallestRealPart(τ), F, 1:1))
    @test real(eigenvalues[smallest_real_part]) ≤ minimum(real.(eigenvalues))

    largest_imaginary_part = first(schur_permutation(LargestImaginaryPart(τ), F, 1:1))
    @test imag(eigenvalues[largest_imaginary_part]) ≥ maximum(imag.(eigenvalues))

    smallest_imaginary_part = first(schur_permutation(SmallestImaginaryPart(τ), F, 1:1))
    @test imag(eigenvalues[smallest_imaginary_part]) ≤ minimum(imag.(eigenvalues))

    # Test whether one can reorder the Schur decomp with multiple eigenvalues
    take = 1 : 3; keep = 4 : 10
    schur_sort!(Near(τ), F, take)
    sorted = F.alpha ./ F.beta
    @test maximum(abs.(sorted[take] .- τ)) ≤ minimum(abs.(sorted[keep] .- τ))

    # Test whether one can reorder the Schur decomp with a single eigenvalue
    schur_sort!(SmallestMagnitude(τ), F, 1)
    sorted = F.alpha ./ F.beta
    @test all(abs(sorted[1]) .≤ abs.(sorted[2 : end]))
end

@testset "Jacobi–Davidson" begin
    function myA!(y, x)
        for i = 1 : length(x)
            @inbounds y[i] = sqrt(i) * x[i]
        end
    end

    function myB!(y, x)
        for i = 1 : length(x)
            @inbounds y[i] = x[i] / sqrt(i)
        end
    end

    struct SuperPreconditioner{numT <: Number}
        target::numT
    end

    function LinearAlgebra.ldiv!(p::SuperPreconditioner{T}, x::AbstractVector{T}) where {T<:Number}
        for i = 1 : length(x)
            @inbounds x[i] = x[i] * sqrt(i) / (i - p.target)
        end
        return x
    end

    function LinearAlgebra.ldiv!(y::AbstractVector{T}, p::SuperPreconditioner{T}, x::AbstractVector{T}) where {T<:Number}
        for i = 1 : length(x)
            @inbounds y[i] = x[i] * sqrt(i) / (i - p.target)
        end
        return y
    end

    n = 10_000
    target = Near(5_000.1 + 0.0im)
    A = LinearMap{Float64}(myA!, n; ismutating = true)
    B = LinearMap{Float64}(myB!, n; ismutating = true)
    P = SuperPreconditioner(target.τ)

    @testset "JDQR" begin
        @testset "Verbosity = $(verbosity)" for verbosity = 1:2
            pschur, residuals = jdqr(A,
                                     solver = BiCGStabl(n, max_mv_products = 10, l = 2),
                                     target = target,
                                     pairs = 5,
                                     tolerance = 1e-9,
                                     subspace_dimensions = 10:20,
                                     max_iter = 1000,
                                     verbosity = verbosity
                                     )
            verbosity == 1 && println("\n\n\n")
            λ = sort(real(pschur.values))

            @test λ ≈ sqrt.(2*(4998:0.5:5000)) rtol=1e-10
        end
    end

    @testset "JDQZ" begin
        @testset "Correction solver = $(label)" for (label,solver) in [("BiCGStabl", BiCGStabl(n, max_mv_products = 10, l = 2)),
                                                                       ("GMRES", GMRES(n))]
            @testset "Verbosity = $(verbosity)" for verbosity = 1:2
                pschur, residuals = jdqz(A, B,
                                         solver = solver,
                                         preconditioner = P,
                                         testspace = Harmonic,
                                         target = target,
                                         pairs = 5,
                                         tolerance = 1e-9,
                                         subspace_dimensions = 10:20,
                                         max_iter = 100,
                                         verbosity = verbosity
                                         )
                verbosity == 1 && println("\n\n\n")
                λ = sort(real(pschur.alphas ./ pschur.betas))

                @test λ ≈ 4998:5002 rtol=1e-10
            end
        end
    end
end

max_stars_repo_name: haampie/jacobi-davidson-julia

max_forks_repo_name: haampie/jacobi-davidson-julia


text:
using Test, Random
import Flux: activations

@testset "basic" begin
  @testset "helpers" begin
    @testset "activations" begin
      dummy_model = Chain(x->x.^2, x->x .- 3, x -> tan.(x))
      x = randn(10)
      @test activations(dummy_model, x)[1] == x.^2
      @test activations(dummy_model, x)[2] == (x.^2 .- 3)
      @test activations(dummy_model, x)[3] == tan.(x.^2 .- 3)

      @test activations(Chain(), x) == ()
      @test activations(Chain(identity, x->:foo), x)[2] == :foo # results include `Any` type
    end
  end

  @testset "Chain" begin
    @test_nowarn Chain(Dense(10, 5, σ), Dense(5, 2))(randn(10))
    @test_throws DimensionMismatch Chain(Dense(10, 5, σ),Dense(2, 1))(randn(10))
    # numeric test should be put into testset of corresponding layer
  end

  @testset "Activations" begin
    c = Chain(Dense(3,5,relu), Dense(5,1,relu))
    X = Float32.([1.0; 1.0; 1.0])
    @test_nowarn gradient(()->Flux.activations(c, X)[2][1], params(c))
  end

  @testset "Dense" begin
    @testset "constructors" begin
      @test size(Dense(10, 100).W) == (100, 10)
      @test Dense(rand(100,10), rand(10)).σ == identity

      @test_throws MethodError Dense(10, 10.5)
      @test_throws MethodError Dense(10, 10.5, tanh)
    end
    @testset "dimensions" begin
      @test  length(Dense(10, 5)(randn(10))) == 5
      @test_throws DimensionMismatch Dense(10, 5)(randn(1))
      @test_throws MethodError Dense(10, 5)(1) # avoid broadcasting
      @test_throws MethodError Dense(10, 5).(randn(10)) # avoid broadcasting
      @test size(Dense(10, 5)(randn(10))) == (5,)
      @test size(Dense(10, 5)(randn(10,2))) == (5,2)
      @test size(Dense(10, 5)(randn(10,2,3))) == (5,2,3)
      @test size(Dense(10, 5)(randn(10,2,3,4))) == (5,2,3,4)
    end
    @testset "zeros" begin
      @test Dense(10, 1, identity, initW = ones, initb = zeros)(ones(10,1)) == 10*ones(1, 1)
      @test Dense(10, 1, identity, initW = ones, initb = zeros)(ones(10,2)) == 10*ones(1, 2)
      @test Dense(10, 2, identity, initW = ones, initb = zeros)(ones(10,1)) == 10*ones(2, 1)
      @test Dense(10, 2, identity, initW = ones, initb = zeros)([ones(10,1) 2*ones(10,1)]) == [10 20; 10 20]
      @test Dense(10, 2, identity, initW = ones, bias = false)([ones(10,1) 2*ones(10,1)]) == [10 20; 10 20]
    end
  end

  @testset "Diagonal" begin
    @test length(Flux.Diagonal(10)(randn(10))) == 10
    @test length(Flux.Diagonal(10)(1)) == 10
    @test length(Flux.Diagonal(10)(randn(1))) == 10
    @test_throws DimensionMismatch Flux.Diagonal(10)(randn(2))

    @test Flux.Diagonal(2)([1 2]) == [1 2; 1 2]
    @test Flux.Diagonal(2)([1,2]) == [1,2]
    @test Flux.Diagonal(2)([1 2; 3 4]) == [1 2; 3 4]

    @test Flux.Diagonal(2)(rand(2,3,4)) |> size == (2, 3, 4)
    @test Flux.Diagonal(2,3)(rand(2,3,4)) |> size == (2, 3, 4)
    @test Flux.Diagonal(2,3,4)(rand(2,3,4)) |> size == (2, 3, 4)
    @test Flux.Diagonal(2,3)(rand(2,1,4)) |> size == (2, 3, 4)
  end

  @testset "Maxout" begin
    # Note that the normal common usage of Maxout is as per the docstring
    # These are abnormal constructors used for testing purposes

    @testset "Constructor" begin
      mo = Maxout(() -> identity, 4)
      input = rand(40)
      @test mo(input) == input
    end

    @testset "simple alternatives" begin
      mo = Maxout((x -> x, x -> 2x, x -> 0.5x))
      input = rand(40)
      @test mo(input) == 2*input
    end

    @testset "complex alternatives" begin
      mo = Maxout((x -> [0.5; 0.1]*x, x -> [0.2; 0.7]*x))
      input = [3.0 2.0]
      target = [0.5, 0.7].*input
      @test mo(input) == target
    end

    @testset "params" begin
      mo = Maxout(()->Dense(32, 64), 4)
      ps = params(mo)
      @test length(ps) == 8  #4 alts, each with weight and bias
    end
  end

  @testset "SkipConnection" begin
    @testset "zero sum" begin
      input = randn(10, 10, 10, 10)
      @test SkipConnection(x -> zeros(size(x)), (a,b) -> a + b)(input) == input
    end

    @testset "concat size" begin
      input = randn(10, 2)
      @test size(SkipConnection(Dense(10,10), (a,b) -> cat(a, b, dims = 2))(input)) == (10,4)
    end
  end

  @testset "Parallel" begin
    @testset "zero sum" begin
      input = randn(10, 10, 10, 10)
      @test Parallel(+, x -> zeros(size(x)), identity)(input) == input
    end

    @testset "concat size" begin
      input = randn(10, 2)
      @test size(Parallel((a, b) -> cat(a, b; dims=2), Dense(10, 10), identity)(input)) == (10, 4)
    end

    @testset "vararg input" begin
      inputs = randn(10), randn(5), randn(4)
      @test size(Parallel(+, Dense(10, 2), Dense(5, 2), Dense(4, 2))(inputs)) == (2,)
    end
  end
end


text:
using GeneralizedChartParsing
using Test

using LinearAlgebra: norm
using Random:        randperm

###################
### Utils Tests ###
###################

@test GeneralizedChartParsing.findallin([1,3,5], [4,3,2,5,1]) == [5,2,4]

############################
### Named Function Tests ###
############################

f = NamedFunction("goo", x -> x + 2)
@test string(f) == "goo"
@test f(3) == 5

##############################
### Construct Test Grammar ###
##############################

categories       = collect(1:10)
start_categories = [1]

double(x)        = (x, x)
countup(x)       = x <= 9 ? (x, x+1) : nothing
terminate(x)     = string(x)

binary_rules     = [double, countup]
terminal_rules   = [terminate]
depcomp          = identity

#############
### Tests ###
#############

g = Grammar(categories, start_categories, binary_rules, terminal_rules, depcomp)
g = add_random_prob_score(g, :count)

terminals = split("1 2 3 4 1 2")
@test s = parse(g, terminals)[1].count == 2

g = add_score(g, :forest, forest_score(g.scores.prob))
s = parse(g, terminals)[1]
t = s.forest |> rand |> collect |> t -> tree_struct(g, t)
ts = ["[1[1[1[1]][2[2[2]][3[3[3]][4[4]]]]][1[1[1]][2[2]]]]",
      "[1[1[1[1[1]][2[2[2]][3[3[3]][4[4]]]]][1[1]]][2[2]]]"]
@test string(t) in ts

g = Grammar(categories, start_categories, binary_rules, terminal_rules, depcomp)
g = add_score(g, :enum_forest, enum_forest_score)
s = parse(g, terminals)[1]
@test string.(tree_structs(g, s.enum_forest)) == ts

s = parse(g, fill("1", 10))[1]
@test length(tree_structs(g, s.enum_forest)) == s.count == GeneralizedChartParsing.number_binary_trees(10)

g = Grammar(categories, start_categories, binary_rules, terminal_rules, depcomp)
g = add_random_prob_score(g, :count)
g = add_expected_count_score(g, :prob)

ec = parse(g, terminals)[1].exp_counts

@test expected_counts_dataframe(g, ec)[:rule] |> unique ==
    [double, countup, terminate]

#################################
### Grammar From String Tests ###
#################################

str = """
    S --> A A | a
    A --> A A | a
    """

g = grammar_from_string(str)
g = add_score(g, :forest, enum_forest_score)
g = add_random_prob_score(g, :count)
g = add_expected_count_score(g, :prob)

rule_strings = ["S --> A A", "S --> a", "A --> A A", "A --> a"]

@test Set(map(string, g.all_rules)) == Set(rule_strings)

p = parse(g, fill("a", 3))["S"]
@test map(string, tree_structs(g, p.forest)) ==
    ["[S[A[a]][A[A[a]][A[a]]]]", "[S[A[A[a]][A[a]]][A[a]]]"]

@test expected_counts_dataframe(g, p.exp_counts)[:rule] .|> string |> Set ==
    Set(["S --> A A", "A --> A A", "A --> a"])

str = """
    S --> A N
    N --> S B | b
    A --> a
    B --> b
    """

g = grammar_from_string(str)
g = add_score(g, :forest, enum_forest_score)

p = parse(g, split("a a b b"))["S"]

@test map(string, tree_structs(g, p.forest)) ==
    ["[S[A[a]][N[S[A[a]][N[b]]][B[b]]]]"]

###############################
### Variational Bayes Tests ###
###############################

g = grammar_from_string("S --> S S | a | b")
g = set_scores(g, (prior = g.scores.count,))
g = add_random_prob_score(g, :prior)
g = add_forest_score(g, :prob)

dataset = let terminals = [fill("a", 3); fill("b", 7)]
    [terminals[randperm(10)] for i in 1:100]
end

loglikelihoods = Float64[]

for i in 1:5
    global g, p = train_grammar(g, dataset)
    push!(loglikelihoods, log(p))
end

true_parameters =
    [9 / 19, 3 / 19, 7 / 19]
inferred_parameters =
    map(float, vec(GeneralizedChartParsing.tabulate_score(g, :prob)))

@test norm(true_parameters - inferred_parameters) < 0.01


text:
export
    laguerre

"""
    laguerre(n::Integer, α::Integer, x::Real)

Generalized Laguerre polynomials

# Examples
```jldoctest
julia> laguerre(0, 1, 5)
1

julia> laguerre(1, 3, 5)
-1
```
"""
function laguerre(n::Integer, α::Integer, x::T) where {T<:Real}
    if n == 0
        return one(T)
    elseif n == 1
        return one(T) + α - x
    else
        L_prev = one(T)
        L_k = one(T) + α - x
        for k = 1:(n-1)
            L_next = ((2k + one(T) + α - x) * L_k - (k+α) * L_prev) / (k+one(T))
            L_prev, L_k = L_k, L_next
        end
        return L_k
    end
end

"""
    laguerre(n::Integer, α::Integer)(x::Real)

Generalized Laguerre polynomials

# Examples
```jldoctest
julia> L = laguerre(1, 3);

julia> L(5)
-1
```
"""
laguerre(n::Integer, α::Integer) = x -> laguerre(n, α, x)

# function hermite(n::T, x::Real) where {T <: Integer}
#     result = 0
#     for k in 0:floor(T, n/2)
#         coeff = (-1) ^ k *
#             2 ^ (n-2k) *
#             factorial(n) / (factorial(k)*factorial(n-2k))

#         x_power = x^(n-2k)

#         result += coeff * x_power
#     end

#     return result
# end


text:
include("./FieldFactory/rayclassgrp.jl")
include("./FieldFactory/ab_exts.jl")
include("./FieldFactory/conductors_generic.jl")
include("./FieldFactory/relative_to_absolute.jl")


text:
#Tree{T} and Node{T} are two parametric abstract types where implementations of the methods could be
#used at multiple places
abstract Tree{T}
abstract Node{T}
abstract Branch{T} <: Node{T} #Branch represents nonterminal nodes
abstract Leaf{T} <: Node{T}

function valid_tree{T}(tree::Tree{T})
   @assert typeof(tree.root)<:Node
   @assert valid_node(tree.root)
   true
end

function valid_node{T}(branch::Branch{T}) #a general definition for noninternal nodes
   @assert valid_node(branch.left)
   @assert valid_node(branch.right)
   @assert typeof(branch.feature)<:Int
   @assert typeof(branch.value)<:Float64
   true
end

function valid_node{T}(leaf::Leaf{T})
   true
end

#calculate the depth of a tree
depth{T}(tree::Tree{T}) = depth(tree.root)
depth{T}(branch::Branch{T}) = 1 + max(depth(branch.left),depth(branch.right))
depth{T}(leaf::Leaf{T}) = 1

depth{T}(tree::Tree{T},node::Node{T}) = depth(tree.root,node)

function depth{T}(branch::Branch{T},node::Node{T})
  if node==branch
    return 1
  end
  left_depth = depth(branch.left,node)
  right_depth = depth(branch.right,node)
  left_depth = left_depth>0?left_depth+1:0 #when left.depth>0, it means that it was in left branch
  right_depth = right_depth>0?right_depth+1:0
  max(left_depth,right_depth)
end
depth{T}(leaf::Leaf{T},node::Node{T}) = node==leaf?1:0

Base.length{T}(tree::Tree{T}) = length(tree.root)
Base.length{T}(branch::Branch{T}) = 1+length(branch.left)+length(branch.right)
Base.length{T}(leaf::Leaf{T}) = 1

parent{T}(tree::Tree{T},node::Node{T}) = parent(tree.root,node)
function parent{T}(branch::Branch{T},node::Node{T})
  a_or_b(a,b) = a==nothing?b:a
  if branch.left==node || branch.right==node
      parentNode = branch
  else
      parentNode = a_or_b(parent(branch.left,node),parent(branch.right,node))
  end
  parentNode
end

parent{T}(leaf::Leaf{T},node::Node{T}) = nothing

function leaves{T}(branch::Branch{T})
  function leaves!{T}(branch::Branch{T},leaf_nodes::Vector{Leaf{T}})
    leaves!(branch.left,leaf_nodes)
    leaves!(branch.right,leaf_nodes)
  end
  function leaves!{T}(leaf::Leaf{T},leaf_nodes::Vector{Leaf{T}})
    push!(leaf_nodes,leaf)
  end
  leaf_nodes = Leaf{T}[]
  leaves!(branch,leaf_nodes)
  leaf_nodes
end

leaves{T}(tree::Tree{T}) = leaves(tree.root)
leaves{T}(leaf::Leaf{T}) = Leaf{T}[leaf];

function branches{T}(branch::Branch{T})
  function branches!{T}(branch::Branch{T}, branch_nodes::Vector{Branch{T}})
    push!(branch_nodes,branch)
    branches!(branch.left,branch_nodes)
    branches!(branch.right,branch_nodes)
  end
  branches!{T}(leaf::Leaf{T},branch_nodes::Vector{Branch{T}}) = nothing
  branch_nodes = Branch{T}[]
  branches!(branch,branch_nodes)
  branch_nodes
end

branches{T}(tree::Tree{T}) = branches(tree.root)
branches{T}(leaf::Leaf{T}) = Branch{T}[]

#grand_branch_nodes are those nodes which have at least one child being a branch
grand_branch{T}(branch::Branch{T}) = typeof(branch.left)<:Branch{T} || typeof(branch.right)<:Branch{T}

function grand_branches{T}(branch::Branch{T})
    function grand_branches!{T}(branch::Branch{T},grand_branch_nodes::Vector{Branch{T}})
      if grand_branch(branch)
        push!(grand_branch_nodes,branch)
        grand_branches!(branch.left,grand_branch_nodes)
        grand_branches!(branch.right,grand_branch_nodes)
      end
    end
    grand_branches!{T}(leaf::Leaf{T},branch_nodes::Vector{Branch{T}}) = nothing

    grand_branch_nodes = Branch{T}[]
    grand_branches!(branch,grand_branch_nodes)
    grand_branch_nodes
end

grand_branches{T}(tree::Tree{T}) = grand_branches(tree.root)
grand_branches{T}(leaf::Leaf{T}) = Branch{T}[]

#not_gran_branch_nodes are those nodes which have exactly two leaves
function not_grand_branches{T}(branch::Branch{T})
  function not_grand_branches!{T}(branch::Branch{T},not_grand_branch_nodes::Vector{Branch{T}})
    if !grand_branch(branch)
      push!(not_grand_branch_nodes,branch)
    else
      not_grand_branches!(branch.left,not_grand_branch_nodes)
      not_grand_branches!(branch.right,not_grand_branch_nodes)
    end
  end

  not_grand_branches!{T}(leaf::Leaf{T},not_grand_branch_nodes::Vector{Branch{T}}) = nothing

  not_grand_branch_nodes = Branch{T}[]
  not_grand_branches!(branch,not_grand_branch_nodes)
  not_grand_branch_nodes
end

not_grand_branches{T}(tree::Tree{T}) = not_grand_branches(tree.root)
not_grand_branches{T}(leaf::Leaf{T}) = Branch{T}[]

macro make_tree_type(Tree,Node)
  eval(quote
       depth(tree::$Tree)              = depth(tree.tree)
       leaves(tree::$Tree)             = leaves(tree.tree)
       branches(tree::$Tree)           = branches(tree.tree)
       grand_branches(tree::$Tree)     = grand_branches(tree.tree)
       not_grand_branches(tree::$Tree) = not_grand_branches(tree.tree)
       depth(tree::$Tree, node::$Node) = depth(tree.tree,node)
       parent(tree::$Tree,node::$Node) = parent(tree.tree,node)
       Base.length(tree::$Tree)        = Base.length(tree.tree)
     end)
end

max_stars_repo_name: mathcg/BayesTree

max_issues_repo_name: mathcg/BayesTree

max_forks_repo_name: mathcg/BayesTree


text:
using CSV, DataFrames, Econometrics, StatsModels, Statistics
sp500 = CSV.read("../Data/sp500.csv", DataFrame)
y = (sp500.rets).^2
n = size(y,1)
x = [ones(n) lags(y,2)]
y = y[3:end]
x = x[3:end,:]
b, varb, e, junk, junk = ols(y,x,vc="nw")
e = [e lags(e,8)]
cor(e[9:end,:])

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
# playLogicalWorld.jl

using Test

##english grammer
include("egrammer.jl")

#Magical world
const QB=:QB
const Homura=:Homura
const Madoka=:Madoka
const Sayaka=:Sayaka
const Kamijo=:Kamijo
const Kyoko=:Kyoko
const Mami=:Mami
const Hitomi=:Hitomi
const Walpurgis=:Walpurgis
const mahosyojo=:mahosyojo
const witch=:witch
const soulgem=:soulgem

MagicalWorld=[
  [Homura love Madoka],
  [Kyoko love Sayaka],
  [Sayaka love Kamijo],
  [Hitomi love Kamijo],
  [Kamijo love Hitomi],
  [Madoka is mahosyojo],
  [Homura is mahosyojo],
  [Kyoko is mahosyojo],
  [Mami is mahosyojo],
  [Sayaka is mahosyojo],
  [Homura have soulgem],
  [Sayaka have soulgem],
  [Kyoko have soulgem],
  [Mami have soulgem],
  [QB make mahosyojo],
  [mahosyojo become witch]
]

MagicalWorld2=[
  [Homura love Madoka],
  [Kyoko love Sayaka],
  [Sayaka love Kamijo],
  [Hitomi love Kamijo],
  [Madoka is mahosyojo],
  [Homura is mahosyojo],
  [Kyoko is mahosyojo],
  [Mami is mahosyojo],
  [Sayaka is mahosyojo],
  [Homura have soulgem],
  [Sayaka have soulgem],
  [Kyoko have soulgem],
  [Mami have soulgem],
  [QB make mahosyojo],
  [mahosyojo become witch]
]

#==
## Query samples
solveQ([who],[[who is mahosyojo]],MagicalWorld)
solveQ([who],[[who have soulgem]],MagicalWorld)
solveQ([who],[[who have soulgem],[who love Kamijo]],MagicalWorld)

solveQ([who,what],[[what become witch]],MagicalWorld)
mm=solveQ([who,what],[[who have soulgem],[what become witch]],MagicalWorld)
statements([who,what],[[who have soulgem],[what become witch]],mm)

##
@testset "reflection rule" begin
 @test solveQ([:who1,:who2],[[:who1 love :who2],[:who2 love :who1]], MagicalWorld)==[[:Hitomi,:Kamijo],[:Kamijo,:Hitomi]]

# [] means no model valid
 @test solveQ([:who1,:who2],[[:who1 love :who2],[:who2 love :who1]], MagicalWorld2)==[]

end
==#

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
function solve()
    counter = Dict("AC"=>0, "WA"=>0, "TLE"=>0, "RE"=>0)
    n = parse(Int, readline())
    for i=1:n
        s = readline()
        counter[s] += 1
    end
    println("AC x $(counter["AC"])")
    println("WA x $(counter["WA"])")
    println("TLE x $(counter["TLE"])")
    println("RE x $(counter["RE"])")
end

solve()

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
module StaggeredFermion_module
    using LinearAlgebra

    import ..Actions:FermiActionParam,FermiActionParam_Wilson,
                FermiActionParam_WilsonClover,FermiActionParam_Staggered
    import ..Gaugefields:GaugeFields,GaugeFields_1d,SU3GaugeFields,SU2GaugeFields,SU3GaugeFields_1d,SU2GaugeFields_1d,
                staggered_phase,SUn,SU2,SU3,SUNGaugeFields,SUNGaugeFields_1d,SU

    #using ..Fermion
    import ..AbstractFermion:FermionFields,
            Wx!,Wdagx!,clear!,substitute_fermion!,Dx!,Dxplus!,
            set_wing_fermi!,fermion_shift!,fermion_shiftB!,add!,WdagWx!,apply_periodicity

    struct StaggeredFermion <: FermionFields
        NC::Int64
        NX::Int64
        NY::Int64
        NZ::Int64
        NT::Int64
        f::Array{ComplexF64,6}

        mass::Float64
        eps::Float64
        Dirac_operator::String
        MaxCGstep::Int64
        BoundaryCondition::Array{Int8,1}
    end

    function StaggeredFermion(NC,NX,NY,NZ,NT,fparam::FermiActionParam,BoundaryCondition)
        mass = fparam.mass
        eps = fparam.eps
        MaxCGstep = fparam.MaxCGstep
        return StaggeredFermion(NC,NX,NY,NZ,NT,mass,eps,MaxCGstep,BoundaryCondition)
    end

    function StaggeredFermion(NC,NX,NY,NZ,NT,mass,eps,MaxCGstep,BoundaryCondition)#r,hop,eps,MaxCGstep)
        Dirac_operator = "Staggered"
        return StaggeredFermion(NC,NX,NY,NZ,NT,zeros(ComplexF64,NC,NX+2,NY+2,NZ+2,NT+2,1),
            mass,eps,Dirac_operator,MaxCGstep,BoundaryCondition)
    end

    function StaggeredFermion(NC,NX,NY,NZ,NT,mass,eps,fermion,MaxCGstep,BoundaryCondition)
        return StaggeredFermion(NC,NX,NY,NZ,NT,zeros(ComplexF64,NC,NX+2,NY+2,NZ+2,NT+2,1),
                    mass,eps,fermion,MaxCGstep,BoundaryCondition)
    end

    function Base.:*(a::StaggeredFermion,b::StaggeredFermion)
        c = 0.0im
        α=1
        for it=1:a.NT
            for iz=1:a.NZ
                for iy=1:a.NY
                    for ix=1:a.NX
                        @simd for ic=1:a.NC
                            c+= conj(a[ic,ix,iy,iz,it,α])*b[ic,ix,iy,iz,it,α]
                        end
                    end
                end
            end
        end

        return c
    end

    function Base.similar(x::StaggeredFermion)
        return StaggeredFermion(x.NC,x.NX,x.NY,x.NZ,x.NT,
                    x.mass,x.eps,x.Dirac_operator,x.MaxCGstep,x.BoundaryCondition)
    end

    function clear!(x::StaggeredFermion,evensite)
        ibush = ifelse(evensite,0,1)
        for it=1:x.NT
            for iz=1:x.NZ
                for iy=1:x.NY
                    xran =1+(1+ibush+iy+iz+it)%2:2:x.NX
                    for ix in xran
                        @simd for ic=1:x.NC
                            x[ic,ix,iy,iz,it,1] = 0
                        end
                    end
                end
            end
        end
        return
    end

    function substitute_fermion!(H,j,x::StaggeredFermion)
        i = 0
        for ialpha = 1:1
            for it=1:x.NT
                for iz=1:x.NZ
                    for iy=1:x.NY
                        for ix=1:x.NX
                            @simd for ic=1:x.NC
                                i += 1
                                H[i,j] = x[ic,ix,iy,iz,it,ialpha]
                            end
                        end
                    end
                end
            end
        end
    end

    function substitute_fermion!(y::Vector,x::StaggeredFermion)
        i = 0
        for ialpha = 1:1
            for it=1:x.NT
                for iz=1:x.NZ
                    for iy=1:x.NY
                        for ix=1:x.NX
                            @simd for ic=1:x.NC
                                i += 1
                                y[i] = x[ic,ix,iy,iz,it,ialpha]
                            end
                        end
                    end
                end
            end
        end
    end

    """
    (D + m)*x
    """
    function Wx!(xout::StaggeredFermion,U::Array{G,1},
        x::StaggeredFermion,temps::Array{T,1},fparam::FermiActionParam_Staggered) where  {T <: FermionFields,G <: GaugeFields}

        temp = temps[4]
        Dx!(temp,U,x,temps,fparam)
        clear!(xout)
        add!(xout,fparam.mass,x,1,temp)
        set_wing_fermi!(xout)
        return

    end

    function Wx!(xout::StaggeredFermion,U::Array{G,1},
        x::StaggeredFermion,temps::Array{T,1},fparam::FermiActionParam_Staggered,indices) where  {T <: FermionFields,G <: GaugeFields}

        temp = temps[4]
        vec_indices1,vec_indices2 = Dx!(temp,U,x,temps,indices)
        clear!(xout)
        add!(xout,fparam.mass,x,1,temp)
        set_wing_fermi!(xout)
        return vec_indices1,vec_indices2

    end

    function Wx!(xout::StaggeredFermion,evensite,U::Array{G,1},
        x::StaggeredFermion,temps::Array{T,1},fparam::FermiActionParam_Staggered) where  {T <: FermionFields,G <: GaugeFields}

        temp = temps[4]
        Dx!(temp,evensite,U,x,temps,fparam)
        clear!(xout)
        add!(xout,fparam.mass,x,1,temp)
        set_wing_fermi!(xout)
        return

    end

    function Dx!(xout::StaggeredFermion,U::Array{G,1},
        x::StaggeredFermion,temps::Array{T,1}) where  {T <: FermionFields,G <: GaugeFields}
        #temp = temps[4]
        temp1 = temps[1]
        temp2 = temps[2]

        #clear!(temp)
        set_wing_fermi!(x)
        clear!(xout)
        for ν=1:4

            fermion_shift!(temp1,U,ν,x)

            fermion_shift!(temp2,U,-ν,x)

            add!(xout,0.5,temp1,-0.5,temp2)

        end
        set_wing_fermi!(xout)

        return
    end

    function Dxplus!(xout::StaggeredFermion,ν::Int64,U::Array{G,1},
        x::StaggeredFermion,temps::Array{T,1}) where  {T <: FermionFields,G <: GaugeFields}
        #temp = temps[4]
        temp1 = temps[1]
        temp2 = temps[2]

        #clear!(temp)
        set_wing_fermi!(x)

        fermion_shift!(temp1,U,ν,x)

        fermion_shift!(temp2,U,-ν,x)

        add!(xout,0.5,temp1,-0.5,temp2)

        set_wing_fermi!(xout)

        return
    end

    function Dx!(xout::StaggeredFermion,U::Array{G,1},
        x::StaggeredFermion,temps::Array{T,1},indices) where  {T <: FermionFields,G <: GaugeFields}
        #temp = temps[4]
        temp1 = temps[1]
        temp2 = temps[2]

        #clear!(temp)
        set_wing_fermi!(x)
        clear!(xout)
        vec_indices1 = Array{Tuple,1}(undef,4)
        vec_indices2 = Array{Tuple,1}(undef,4)

        for ν=1:4
            clear!(temp1)
            vec_indices1[ν] = fermion_shift!(temp1,U,ν,x,indices)

            #println(temp1*temp1)

            clear!(temp2)

            vec_indices2[ν] = fermion_shift!(temp2,U,-ν,x,indices)

            add!(xout,0.5,temp1,-0.5,temp2)

        end
        set_wing_fermi!(xout)

        return vec_indices1,vec_indices2
    end

    function DdagDx!(xout::StaggeredFermion,U::Array{G,1},
        x::StaggeredFermion,temps::Array{T,1},indices) where  {T <: FermionFields,G <: GaugeFields}
        #temp = temps[4]

        temp1 = temps[1]
        temp2 = temps[2]
        temp3 = temps[3]
        temp4 = temps[4]
        temp5 = temps[5]

        clear!(temp1)
        clear!(temp2)
        clear!(temp4)
        clear!(temp5)

        #clear!(temp)
        set_wing_fermi!(x)
        clear!(xout)
        vec_indices = Array{Tuple,1}(undef,8)
        #vec_indices2 = Array{Tuple,1}(undef,4)

        clear!(temp3)
        for ν=1:4
            vec_indices[2ν-1] = fermion_shift!(temp1,U,ν,x,[indices])[1]
            vec_indices[2ν] = fermion_shift!(temp2,U,-ν,x,[indices])[1]

        end
        add!(temp3,0.5,temp1,-0.5,temp2)

        for ν=1:4
            fermion_shift!(temp4,U,ν,temp3,vec_indices)
            fermion_shift!(temp5,U,-ν,temp3,vec_indices)

        end
        add!(xout,0.5,temp4,-0.5,temp5)

        set_wing_fermi!(xout)
        return

    end

    function Dx!(xout::StaggeredFermion,U::Array{G,1},
        x::StaggeredFermion,temps::Array{T,1},fparam::FermiActionParam_Staggered) where  {T <: FermionFields,G <: GaugeFields}
        Dx!(xout,U,x,temps)
        return
    end

    function Dx!(xout::StaggeredFermion,evensite,U::Array{G,1},
        x::StaggeredFermion,temps::Array{T,1},fparam::FermiActionParam_Staggered) where  {T <: FermionFields,G <: GaugeFields}
        #temp = temps[4]
        temp1 = temps[1]
        temp2 = temps[2]

        #clear!(temp)
        set_wing_fermi!(x)
        clear!(xout)
        for ν=1:4
            clear!(temp1)
            fermion_shift!(temp1,evensite,U,ν,x)

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

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

            #add!(temp,1,temp1,-1,temp2)
            add!(xout,0.5,temp1,-0.5,temp2)

        end
        set_wing_fermi!(xout)

        #clear!(xout)
        #add!(xout,fparam.mass,x,1,temp)

        #display(xout)
        #    exit()
        return
    end

    function Wdagx!(xout::StaggeredFermion,U::Array{G,1},
        x::StaggeredFermion,temps::Array{T,1},fparam::FermiActionParam_Staggered) where {T <: FermionFields,G <: GaugeFields}

        temp = temps[4]
        Dx!(temp,U,x,temps,fparam)
        clear!(xout)
        add!(xout,fparam.mass,x,-1,temp)
        return

    end

    function Wdagx!(xout::StaggeredFermion,U::Array{G,1},
        x::StaggeredFermion,temps::Array{T,1},fparam::FermiActionParam_Staggered,vec_indices1,vec_indices2) where {T <: FermionFields,G <: GaugeFields}

        temp = temps[4]
        Dx!(temp,U,x,temps,vec_indices1,vec_indices2)
        clear!(xout)
        add!(xout,fparam.mass,x,-1,temp)
        return

    end

    function WdagWx!(xout::T,U::Array{G,1},x::T,temps::Array{T,1},mass::Number,indices) where {T <: StaggeredFermion,G <: GaugeFields}

        temp = temps[6]
        #println("x = ", x*x)
        DdagDx!(temp,U,x,temps,indices)

        #println("temp ",temp*temp)
        clear!(xout)
        add!(xout,mass^2,x,-1,temp)
        #println("xout ",xout*xout)
        return
    end

    function WdagWx!(xout::T,U::Array{G,1},x::T,temps::Array{T,1},fparam::FermiActionParam_Staggered,indices) where {T <: StaggeredFermion,G <: GaugeFields}
        WdagWx!(xout,U,x,temps,fparam.mass,indices)
        return
    end

    function Ddagx!(xout::StaggeredFermion,U::Array{G,1},
        x::StaggeredFermion,temps::Array{T,1},fparam::FermiActionParam_Staggered) where  {T <: FermionFields,G <: GaugeFields}
        #temp = temps[4]
        temp1 = temps[1]
        temp2 = temps[2]

        #clear!(temp)
        clear!(xout)
        set_wing_fermi!(x)
        for ν=1:4
            fermion_shift!(temp1,U,ν,x)

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

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

            #add!(temp,1,temp1,-1,temp2)

            add!(xout,-0.5,temp1,0.5,temp2)
            #add!(xout,1,temp1,-1,temp2)

        end
        set_wing_fermi!(xout)

        #clear!(xout)
        #add!(xout,fparam.mass,x,1,temp)

        #display(xout)
        #    exit()
        return
    end

    function set_wing_fermi!(a::StaggeredFermion)
        NT = a.NT
        NZ = a.NZ
        NY = a.NY
        NX = a.NX
        NC = a.NC

        #!  X-direction
        for ialpha=1:1
            for it=1:NT
                for iz = 1:NZ
                    for iy=1:NY
                        @simd for k=1:NC
                            a[k,0,iy,iz,it,ialpha] = a.BoundaryCondition[1]*a[k,NX,iy,iz,it,ialpha]
                        end
                    end
                end
            end
        end

        for ialpha=1:1
            for it=1:NT
                for iz=1:NZ
                    for iy=1:NY
                        @simd for k=1:NC
                            a[k,NX+1,iy,iz,it,ialpha] = a.BoundaryCondition[1]*a[k,1,iy,iz,it,ialpha]
                        end
                    end
                end
            end
        end

        #Y-direction
        for ialpha = 1:1
            for it=1:NT
                for iz=1:NZ
                    for ix=1:NX
                        @simd for k=1:NC
                            a[k,ix,0,iz,it,ialpha] = a.BoundaryCondition[2]*a[k,ix,NY,iz,it,ialpha]
                        end
                    end
                end
            end
        end

        for ialpha=1:1
            for it=1:NT
                for iz=1:NZ
                    for ix=1:NX
                        @simd for k=1:NC
                            a[k,ix,NY+1,iz,it,ialpha] = a.BoundaryCondition[2]*a[k,ix,1,iz,it,ialpha]
                        end
                    end
                end
            end
        end

        for ialpha=1:1
            # Z-direction
            for it=1:NT
                for iy=1:NY
                    for ix=1:NX
                        @simd for k=1:NC
                            a[k,ix,iy,0,it,ialpha] = a.BoundaryCondition[3]*a[k,ix,iy,NZ,it,ialpha]
                            a[k,ix,iy,NZ+1,it,ialpha] = a.BoundaryCondition[3]*a[k,ix,iy,1,it,ialpha]

                        end
                    end
                end
            end

            #T-direction
            for iz=1:NZ
                for iy=1:NY
                    for ix=1:NX
                        @simd for k=1:NC
                            a[k,ix,iy,iz,0,ialpha] = a.BoundaryCondition[4]*a[k,ix,iy,iz,NT,ialpha]
                            a[k,ix,iy,iz,NT+1,ialpha] = a.BoundaryCondition[4]*a[k,ix,iy,iz,1,ialpha]
                        end
                    end
                end
            end

        end

    end

    function fermion_shift!(b::StaggeredFermion,u::Array{T,1},μ::Int,a::StaggeredFermion) where T <: SU3GaugeFields
        if μ == 0
            substitute!(b,a)
            return
        end

        NX = a.NX
        NY = a.NY
        NZ = a.NZ
        NT = a.NT
        NC = a.NC

        #NTrange = get_looprange(NT)
        #println(NTrange)
        if μ > 0
            #idel = zeros(Int64,4)
            #idel[μ] = 1

            n6 = size(a.f)[6]
            for ialpha=1:1
                #for it=NTrange
                for it=1:NT
                    it1 = it + ifelse(μ ==4,1,0) #idel[4]
                    for iz=1:NZ
                        iz1 = iz + ifelse(μ ==3,1,0) #idel[3]
                        for iy=1:NY
                            iy1 = iy + ifelse(μ ==2,1,0) #idel[2]
                            @simd for ix=1:NX
                                ix1 = ix + ifelse(μ ==1,1,0) #idel[1]
                                η = staggered_phase(μ,ix,iy,iz,it,NX,NY,NZ,NT)

                                b[1,ix,iy,iz,it,ialpha] = η*(u[μ][1,1,ix,iy,iz,it]*a[1,ix1,iy1,iz1,it1,ialpha] +
                                                            u[μ][1,2,ix,iy,iz,it]*a[2,ix1,iy1,iz1,it1,ialpha] +
                                                            u[μ][1,3,ix,iy,iz,it]*a[3,ix1,iy1,iz1,it1,ialpha])

                                b[2,ix,iy,iz,it,ialpha] = η*(u[μ][2,1,ix,iy,iz,it]*a[1,ix1,iy1,iz1,it1,ialpha] +
                                                            u[μ][2,2,ix,iy,iz,it]*a[2,ix1,iy1,iz1,it1,ialpha] +
                                                            u[μ][2,3,ix,iy,iz,it]*a[3,ix1,iy1,iz1,it1,ialpha])

                                b[3,ix,iy,iz,it,ialpha] = η*(u[μ][3,1,ix,iy,iz,it]*a[1,ix1,iy1,iz1,it1,ialpha] +
                                                            u[μ][3,2,ix,iy,iz,it]*a[2,ix1,iy1,iz1,it1,ialpha] +
                                                            u[μ][3,3,ix,iy,iz,it]*a[3,ix1,iy1,iz1,it1,ialpha])

                                #=
                                for k1=1:3
                                    b[k1,ix,iy,iz,it,ialpha] = 0
                                    for k2=1:3
                                        b[k1,ix,iy,iz,it,ialpha] += u[μ][k1,k2,ix,iy,iz,it]*a[k2,ix1,iy1,iz1,it1,ialpha]
                                    end
                                end
                                =#
                            end
                        end
                    end
                end
            end

        elseif μ < 0
            #idel = zeros(Int64,4)
            #idel[-μ] = 1
            #n6 = size(b.f)[6]
            for ialpha =1:1
                for it=1:NT
                    it1 = it - ifelse(-μ ==4,1,0) #idel[4]
                    for iz=1:NZ
                        iz1 = iz - ifelse(-μ ==3,1,0) #idel[3]
                        for iy=1:NY
                            iy1 = iy - ifelse(-μ ==2,1,0)  #idel[2]
                            @simd for ix=1:NX
                                ix1 = ix - ifelse(-μ ==1,1,0) #idel[1]

                                η = staggered_phase(-μ,ix1,iy1,iz1,it1,NX,NY,NZ,NT)

                                b[1,ix,iy,iz,it,ialpha] = η*(conj(u[-μ][1,1,ix1,iy1,iz1,it1])*a[1,ix1,iy1,iz1,it1,ialpha] +
                                                            conj(u[-μ][2,1,ix1,iy1,iz1,it1])*a[2,ix1,iy1,iz1,it1,ialpha] +
                                                            conj(u[-μ][3,1,ix1,iy1,iz1,it1])*a[3,ix1,iy1,iz1,it1,ialpha])

                                b[2,ix,iy,iz,it,ialpha] = η*(conj(u[-μ][1,2,ix1,iy1,iz1,it1])*a[1,ix1,iy1,iz1,it1,ialpha] +
                                                            conj(u[-μ][2,2,ix1,iy1,iz1,it1])*a[2,ix1,iy1,iz1,it1,ialpha] +
                                                            conj(u[-μ][3,2,ix1,iy1,iz1,it1])*a[3,ix1,iy1,iz1,it1,ialpha])

                                b[3,ix,iy,iz,it,ialpha] = η*(conj(u[-μ][1,3,ix1,iy1,iz1,it1])*a[1,ix1,iy1,iz1,it1,ialpha] +
                                                            conj(u[-μ][2,3,ix1,iy1,iz1,it1])*a[2,ix1,iy1,iz1,it1,ialpha] +
                                                            conj(u[-μ][3,3,ix1,iy1,iz1,it1])*a[3,ix1,iy1,iz1,it1,ialpha])
                                #=
                                for k1=1:3
                                    b[k1,ix,iy,iz,it,ialpha] = 0
                                    for k2=1:3
                                        b[k1,ix,iy,iz,it,ialpha] += conj(u[-μ][k2,k1,ix1,iy1,iz1,it1])*a[k2,ix1,iy1,iz1,it1,ialpha]
                                    end
                                end
                                =#
                            end
                        end
                    end
                end
            end
        end

    end

    function fermion_shift!(b::StaggeredFermion,u::Array{GaugeFields{SU{NC}},1},μ::Int,a::StaggeredFermion) where NC
        if μ == 0
            substitute!(b,a)
            return
        end

        NX = a.NX
        NY = a.NY
        NZ = a.NZ
        NT = a.NT
        #NC = a.NC

        #NTrange = get_looprange(NT)
        #println(NTrange)
        if μ > 0
            #idel = zeros(Int64,4)
            #idel[μ] = 1

            n6 = size(a.f)[6]
            for ialpha=1:1
                #for it=NTrange
                for it=1:NT
                    it1 = it + ifelse(μ ==4,1,0) #idel[4]
                    for iz=1:NZ
                        iz1 = iz + ifelse(μ ==3,1,0) #idel[3]
                        for iy=1:NY
                            iy1 = iy + ifelse(μ ==2,1,0) #idel[2]
                            for ix=1:NX
                                ix1 = ix + ifelse(μ ==1,1,0) #idel[1]
                                η = staggered_phase(μ,ix,iy,iz,it,NX,NY,NZ,NT)

                                for k1=1:NC
                                    b[k1,ix,iy,iz,it,ialpha] = 0
                                    @simd for k2=1:NC
                                        b[k1,ix,iy,iz,it,ialpha] += η*u[μ][k1,k2,ix,iy,iz,it]*a[k2,ix1,iy1,iz1,it1,ialpha]
                                    end
                                end

                            end
                        end
                    end
                end
            end

        elseif μ < 0
            #idel = zeros(Int64,4)
            #idel[-μ] = 1
            #n6 = size(b.f)[6]
            for ialpha =1:1
                for it=1:NT
                    it1 = it - ifelse(-μ ==4,1,0) #idel[4]
                    for iz=1:NZ
                        iz1 = iz - ifelse(-μ ==3,1,0) #idel[3]
                        for iy=1:NY
                            iy1 = iy - ifelse(-μ ==2,1,0)  #idel[2]
                            for ix=1:NX
                                ix1 = ix - ifelse(-μ ==1,1,0) #idel[1]

                                η = staggered_phase(-μ,ix1,iy1,iz1,it1,NX,NY,NZ,NT)

                                for k1=1:NC
                                    b[k1,ix,iy,iz,it,ialpha] = 0
                                    @simd for k2=1:NC
                                        b[k1,ix,iy,iz,it,ialpha] += η*conj(u[-μ][k2,k1,ix1,iy1,iz1,it1])*a[k2,ix1,iy1,iz1,it1,ialpha]
                                    end
                                end

                            end
                        end
                    end
                end
            end
        end

    end

    function fermion_shift!(b::StaggeredFermion,evensite,u::Array{T,1},μ::Int,a::StaggeredFermion) where T <: SU3GaugeFields
        if μ == 0
            substitute!(b,a)
            return
        end

        ibush = ifelse(evensite,0,1)

        NX = a.NX
        NY = a.NY
        NZ = a.NZ
        NT = a.NT
        NC = a.NC

        #NTrange = get_looprange(NT)
        #println(NTrange)
        if μ > 0
            #idel = zeros(Int64,4)
            #idel[μ] = 1

            n6 = size(a.f)[6]
            for ialpha=1:1
                #for it=NTrange
                for it=1:NT
                    it1 = it + ifelse(μ ==4,1,0) #idel[4]
                    for iz=1:NZ
                        iz1 = iz + ifelse(μ ==3,1,0) #idel[3]
                        for iy=1:NY
                            iy1 = iy + ifelse(μ ==2,1,0) #idel[2]
                            xran =1+(1+ibush+iy+iz+it)%2:2:NX
                            @simd for ix in xran
                            #for ix=1:NX
                                ix1 = ix + ifelse(μ ==1,1,0) #idel[1]
                                η = staggered_phase(μ,ix,iy,iz,it,NX,NY,NZ,NT)

                                b[1,ix,iy,iz,it,ialpha] = η*(u[μ][1,1,ix,iy,iz,it]*a[1,ix1,iy1,iz1,it1,ialpha] +
                                                            u[μ][1,2,ix,iy,iz,it]*a[2,ix1,iy1,iz1,it1,ialpha] +
                                                            u[μ][1,3,ix,iy,iz,it]*a[3,ix1,iy1,iz1,it1,ialpha])

                                b[2,ix,iy,iz,it,ialpha] = η*(u[μ][2,1,ix,iy,iz,it]*a[1,ix1,iy1,iz1,it1,ialpha] +
                                                            u[μ][2,2,ix,iy,iz,it]*a[2,ix1,iy1,iz1,it1,ialpha] +
                                                            u[μ][2,3,ix,iy,iz,it]*a[3,ix1,iy1,iz1,it1,ialpha])

                                b[3,ix,iy,iz,it,ialpha] = η*(u[μ][3,1,ix,iy,iz,it]*a[1,ix1,iy1,iz1,it1,ialpha] +
                                                            u[μ][3,2,ix,iy,iz,it]*a[2,ix1,iy1,iz1,it1,ialpha] +
                                                            u[μ][3,3,ix,iy,iz,it]*a[3,ix1,iy1,iz1,it1,ialpha])

                                #=
                                for k1=1:3
                                    b[k1,ix,iy,iz,it,ialpha] = 0
                                    for k2=1:3
                                        b[k1,ix,iy,iz,it,ialpha] += u[μ][k1,k2,ix,iy,iz,it]*a[k2,ix1,iy1,iz1,it1,ialpha]
                                    end
                                end
                                =#
                            end
                        end
                    end
                end
            end

        elseif μ < 0
            #idel = zeros(Int64,4)
            #idel[-μ] = 1
            #n6 = size(b.f)[6]
            for ialpha =1:1
                for it=1:NT
                    it1 = it - ifelse(-μ ==4,1,0) #idel[4]
                    for iz=1:NZ
                        iz1 = iz - ifelse(-μ ==3,1,0) #idel[3]
                        for iy=1:NY
                            iy1 = iy - ifelse(-μ ==2,1,0)  #idel[2]
                            xran =1+(1+ibush+iy+iz+it)%2:2:NX
                            @simd for ix in xran
                                ix1 = ix - ifelse(-μ ==1,1,0) #idel[1]

                                η = staggered_phase(-μ,ix1,iy1,iz1,it1,NX,NY,NZ,NT)

                                b[1,ix,iy,iz,it,ialpha] = η*(conj(u[-μ][1,1,ix1,iy1,iz1,it1])*a[1,ix1,iy1,iz1,it1,ialpha] +
                                                            conj(u[-μ][2,1,ix1,iy1,iz1,it1])*a[2,ix1,iy1,iz1,it1,ialpha] +
                                                            conj(u[-μ][3,1,ix1,iy1,iz1,it1])*a[3,ix1,iy1,iz1,it1,ialpha])

                                b[2,ix,iy,iz,it,ialpha] = η*(conj(u[-μ][1,2,ix1,iy1,iz1,it1])*a[1,ix1,iy1,iz1,it1,ialpha] +
                                                            conj(u[-μ][2,2,ix1,iy1,iz1,it1])*a[2,ix1,iy1,iz1,it1,ialpha] +
                                                            conj(u[-μ][3,2,ix1,iy1,iz1,it1])*a[3,ix1,iy1,iz1,it1,ialpha])

                                b[3,ix,iy,iz,it,ialpha] = η*(conj(u[-μ][1,3,ix1,iy1,iz1,it1])*a[1,ix1,iy1,iz1,it1,ialpha] +
                                                            conj(u[-μ][2,3,ix1,iy1,iz1,it1])*a[2,ix1,iy1,iz1,it1,ialpha] +
                                                            conj(u[-μ][3,3,ix1,iy1,iz1,it1])*a[3,ix1,iy1,iz1,it1,ialpha])
                                #=
                                for k1=1:3
                                    b[k1,ix,iy,iz,it,ialpha] = 0
                                    for k2=1:3
                                        b[k1,ix,iy,iz,it,ialpha] += conj(u[-μ][k2,k1,ix1,iy1,iz1,it1])*a[k2,ix1,iy1,iz1,it1,ialpha]
                                    end
                                end
                                =#
                            end
                        end
                    end
                end
            end
        end

    end

    function fermion_shift!(b::StaggeredFermion,u::Array{T,1},μ::Int,a::StaggeredFermion) where T <: SU2GaugeFields
        if μ == 0
            substitute!(b,a)
            return
        end

        NX = a.NX
        NY = a.NY
        NZ = a.NZ
        NT = a.NT
        NC = a.NC

        #NTrange = get_looprange(NT)
        #println(NTrange)
        if μ > 0
            #idel = zeros(Int64,4)
            #idel[μ] = 1

            n6 = size(a.f)[6]
            for ialpha=1:1
                #for it=NTrange
                for it=1:NT
                    it1 = it + ifelse(μ ==4,1,0) #idel[4]
                    for iz=1:NZ
                        iz1 = iz + ifelse(μ ==3,1,0) #idel[3]
                        for iy=1:NY
                            iy1 = iy + ifelse(μ ==2,1,0) #idel[2]
                            for ix=1:NX
                                ix1 = ix + ifelse(μ ==1,1,0) #idel[1]
                                η = staggered_phase(μ,ix,iy,iz,it,NX,NY,NZ,NT)

                                b[1,ix,iy,iz,it,ialpha] = η*(u[μ][1,1,ix,iy,iz,it]*a[1,ix1,iy1,iz1,it1,ialpha] +
                                                            u[μ][1,2,ix,iy,iz,it]*a[2,ix1,iy1,iz1,it1,ialpha] )

                                b[2,ix,iy,iz,it,ialpha] = η*(u[μ][2,1,ix,iy,iz,it]*a[1,ix1,iy1,iz1,it1,ialpha] +
                                                            u[μ][2,2,ix,iy,iz,it]*a[2,ix1,iy1,iz1,it1,ialpha] )

                                #=
                                for k1=1:3
                                    b[k1,ix,iy,iz,it,ialpha] = 0
                                    for k2=1:3
                                        b[k1,ix,iy,iz,it,ialpha] += u[μ][k1,k2,ix,iy,iz,it]*a[k2,ix1,iy1,iz1,it1,ialpha]
                                    end
                                end
                                =#
                            end
                        end
                    end
                end
            end

        elseif μ < 0
            #idel = zeros(Int64,4)
            #idel[-μ] = 1
            #n6 = size(b.f)[6]
            for ialpha =1:1
                for it=1:NT
                    it1 = it - ifelse(-μ ==4,1,0) #idel[4]
                    for iz=1:NZ
                        iz1 = iz - ifelse(-μ ==3,1,0) #idel[3]
                        for iy=1:NY
                            iy1 = iy - ifelse(-μ ==2,1,0)  #idel[2]
                            for ix=1:NX
                                ix1 = ix - ifelse(-μ ==1,1,0) #idel[1]
                                η = staggered_phase(-μ,ix1,iy1,iz1,it1,NX,NY,NZ,NT)

                                b[1,ix,iy,iz,it,ialpha] = η*(conj(u[-μ][1,1,ix1,iy1,iz1,it1])*a[1,ix1,iy1,iz1,it1,ialpha] +
                                                            conj(u[-μ][2,1,ix1,iy1,iz1,it1])*a[2,ix1,iy1,iz1,it1,ialpha] )

                                b[2,ix,iy,iz,it,ialpha] = η*(conj(u[-μ][1,2,ix1,iy1,iz1,it1])*a[1,ix1,iy1,iz1,it1,ialpha] +
                                                            conj(u[-μ][2,2,ix1,iy1,iz1,it1])*a[2,ix1,iy1,iz1,it1,ialpha] )

                            end
                        end
                    end
                end
            end
        end

    end

    function fermion_shift!(b::StaggeredFermion,u::Array{T,1},μ::Int,a::StaggeredFermion,vec_indices) where T <: SU2GaugeFields
        if μ == 0
            substitute!(b,a)
            return
        end

        NX = a.NX
        NY = a.NY
        NZ = a.NZ
        NT = a.NT
        NC = a.NC

        #NTrange = get_looprange(NT)
        #println(NTrange)
        outindices = Array{Tuple,1}(undef,length(vec_indices))
        #clear!(b)
        if μ > 0
            for (k,indices) in enumerate(vec_indices)

                ix1 = indices[1]
                iy1 = indices[2]
                iz1 = indices[3]
                it1 = indices[4]
                ialpha = indices[5]

                ix = ix1 - ifelse(μ ==1,1,0)
                iy = iy1 - ifelse(μ ==2,1,0)
                iz = iz1 - ifelse(μ ==3,1,0)
                it = it1 - ifelse(μ ==4,1,0)

                η = staggered_phase(μ,ix,iy,iz,it,NX,NY,NZ,NT)

                ix,iy,iz,it,sign =  apply_periodicity(ix,iy,iz,it,NX,NY,NZ,NT,b.BoundaryCondition)

                outindices[k] = (ix,iy,iz,it,ialpha)

                b[1,ix,iy,iz,it,ialpha] += sign*η*(u[μ][1,1,ix,iy,iz,it]*a[1,ix1,iy1,iz1,it1,ialpha] +
                                            u[μ][1,2,ix,iy,iz,it]*a[2,ix1,iy1,iz1,it1,ialpha] )

                b[2,ix,iy,iz,it,ialpha] += sign*η*(u[μ][2,1,ix,iy,iz,it]*a[1,ix1,iy1,iz1,it1,ialpha] +
                                            u[μ][2,2,ix,iy,iz,it]*a[2,ix1,iy1,iz1,it1,ialpha] )
            end

        elseif μ < 0
            for (k,indices) in enumerate(vec_indices)

                ix1 = indices[1]
                iy1 = indices[2]
                iz1 = indices[3]
                it1 = indices[4]
                ialpha = indices[5]

                ix = ix1 + ifelse(-μ ==1,1,0)
                iy = iy1 + ifelse(-μ ==2,1,0)
                iz = iz1 + ifelse(-μ ==3,1,0)
                it = it1 + ifelse(-μ ==4,1,0)

                η = staggered_phase(-μ,ix1,iy1,iz1,it1,NX,NY,NZ,NT)

                ix,iy,iz,it,sign =  apply_periodicity(ix,iy,iz,it,NX,NY,NZ,NT,b.BoundaryCondition)

                outindices[k] = (ix,iy,iz,it,ialpha)

                b[1,ix,iy,iz,it,ialpha] += sign*η*(conj(u[-μ][1,1,ix1,iy1,iz1,it1])*a[1,ix1,iy1,iz1,it1,ialpha] +
                                            conj(u[-μ][2,1,ix1,iy1,iz1,it1])*a[2,ix1,iy1,iz1,it1,ialpha] )

                b[2,ix,iy,iz,it,ialpha] += sign*η*(conj(u[-μ][1,2,ix1,iy1,iz1,it1])*a[1,ix1,iy1,iz1,it1,ialpha] +
                                            conj(u[-μ][2,2,ix1,iy1,iz1,it1])*a[2,ix1,iy1,iz1,it1,ialpha] )
            end

        end

        return outindices

    end

    function fermion_shift!(b::StaggeredFermion,u::Array{T,1},μ::Int,a::StaggeredFermion,vec_indices) where T <: SU3GaugeFields
        if μ == 0
            substitute!(b,a)
            return
        end

        NX = a.NX
        NY = a.NY
        NZ = a.NZ
        NT = a.NT
        NC = a.NC

        #NTrange = get_looprange(NT)
        #println(NTrange)
        outindices = Array{Tuple,1}(undef,length(vec_indices))
        #clear!(b)
        if μ > 0
            for (k,indices) in enumerate(vec_indices)

                ix1 = indices[1]
                iy1 = indices[2]
                iz1 = indices[3]
                it1 = indices[4]
                ialpha = indices[5]

                ix = ix1 - ifelse(μ ==1,1,0)
                iy = iy1 - ifelse(μ ==2,1,0)
                iz = iz1 - ifelse(μ ==3,1,0)
                it = it1 - ifelse(μ ==4,1,0)

                η = staggered_phase(μ,ix,iy,iz,it,NX,NY,NZ,NT)

                ix,iy,iz,it,sign =  apply_periodicity(ix,iy,iz,it,NX,NY,NZ,NT,b.BoundaryCondition)

                outindices[k] = (ix,iy,iz,it,ialpha)

                b[1,ix,iy,iz,it,ialpha] += sign*η*(u[μ][1,1,ix,iy,iz,it]*a[1,ix1,iy1,iz1,it1,ialpha] +
                                            u[μ][1,2,ix,iy,iz,it]*a[2,ix1,iy1,iz1,it1,ialpha] +
                                            u[μ][1,3,ix,iy,iz,it]*a[3,ix1,iy1,iz1,it1,ialpha]
                                            )

                b[2,ix,iy,iz,it,ialpha] += sign*η*(u[μ][2,1,ix,iy,iz,it]*a[1,ix1,iy1,iz1,it1,ialpha] +
                                            u[μ][2,2,ix,iy,iz,it]*a[2,ix1,iy1,iz1,it1,ialpha] +
                                            u[μ][2,3,ix,iy,iz,it]*a[3,ix1,iy1,iz1,it1,ialpha] )
                b[3,ix,iy,iz,it,ialpha] += sign*η*(u[μ][3,1,ix,iy,iz,it]*a[1,ix1,iy1,iz1,it1,ialpha] +
                                            u[μ][3,2,ix,iy,iz,it]*a[2,ix1,iy1,iz1,it1,ialpha] +
                                            u[μ][3,3,ix,iy,iz,it]*a[3,ix1,iy1,iz1,it1,ialpha] )
            end

        elseif μ < 0
            for (k,indices) in enumerate(vec_indices)

                ix1 = indices[1]
                iy1 = indices[2]
                iz1 = indices[3]
                it1 = indices[4]
                ialpha = indices[5]

                ix = ix1 + ifelse(-μ ==1,1,0)
                iy = iy1 + ifelse(-μ ==2,1,0)
                iz = iz1 + ifelse(-μ ==3,1,0)
                it = it1 + ifelse(-μ ==4,1,0)

                η = staggered_phase(-μ,ix1,iy1,iz1,it1,NX,NY,NZ,NT)

                ix,iy,iz,it,sign =  apply_periodicity(ix,iy,iz,it,NX,NY,NZ,NT,b.BoundaryCondition)

                outindices[k] = (ix,iy,iz,it,ialpha)

                b[1,ix,iy,iz,it,ialpha] += sign*η*(conj(u[-μ][1,1,ix1,iy1,iz1,it1])*a[1,ix1,iy1,iz1,it1,ialpha] +
                                            conj(u[-μ][2,1,ix1,iy1,iz1,it1])*a[2,ix1,iy1,iz1,it1,ialpha] +
                                            conj(u[-μ][3,1,ix1,iy1,iz1,it1])*a[3,ix1,iy1,iz1,it1,ialpha] )

                b[2,ix,iy,iz,it,ialpha] += sign*η*(conj(u[-μ][1,2,ix1,iy1,iz1,it1])*a[1,ix1,iy1,iz1,it1,ialpha] +
                                            conj(u[-μ][2,2,ix1,iy1,iz1,it1])*a[2,ix1,iy1,iz1,it1,ialpha]+
                                            conj(u[-μ][3,2,ix1,iy1,iz1,it1])*a[3,ix1,iy1,iz1,it1,ialpha] )
                b[3,ix,iy,iz,it,ialpha] += sign*η*(conj(u[-μ][1,3,ix1,iy1,iz1,it1])*a[1,ix1,iy1,iz1,it1,ialpha] +
                                            conj(u[-μ][2,3,ix1,iy1,iz1,it1])*a[2,ix1,iy1,iz1,it1,ialpha]+
                                            conj(u[-μ][3,3,ix1,iy1,iz1,it1])*a[3,ix1,iy1,iz1,it1,ialpha] )
            end

        end

        return outindices

    end

    function fermion_shift!(b::StaggeredFermion,u::Array{GaugeFields{SU{NC}},1},μ::Int,a::StaggeredFermion,vec_indices) where NC
        if μ == 0
            substitute!(b,a)
            return
        end

        NX = a.NX
        NY = a.NY
        NZ = a.NZ
        NT = a.NT
        #NC = a.NC

        #NTrange = get_looprange(NT)
        #println(NTrange)
        outindices = Array{Tuple,1}(undef,length(vec_indices))
        #clear!(b)
        if μ > 0
            for (k,indices) in enumerate(vec_indices)

                ix1 = indices[1]
                iy1 = indices[2]
                iz1 = indices[3]
                it1 = indices[4]
                ialpha = indices[5]

                ix = ix1 - ifelse(μ ==1,1,0)
                iy = iy1 - ifelse(μ ==2,1,0)
                iz = iz1 - ifelse(μ ==3,1,0)
                it = it1 - ifelse(μ ==4,1,0)

                η = staggered_phase(μ,ix,iy,iz,it,NX,NY,NZ,NT)

                ix,iy,iz,it,sign =  apply_periodicity(ix,iy,iz,it,NX,NY,NZ,NT,b.BoundaryCondition)

                outindices[k] = (ix,iy,iz,it,ialpha)

                for k1=1:NC
                    b[k1,ix,iy,iz,it,ialpha] = 0
                    for k2=1:NC
                        b[k1,ix,iy,iz,it,ialpha] += sign*η*u[μ][k1,k2,ix,iy,iz,it]*a[k2,ix1,iy1,iz1,it1,ialpha]
                    end
                end

            end

        elseif μ < 0
            for (k,indices) in enumerate(vec_indices)

                ix1 = indices[1]
                iy1 = indices[2]
                iz1 = indices[3]
                it1 = indices[4]
                ialpha = indices[5]

                ix = ix1 + ifelse(-μ ==1,1,0)
                iy = iy1 + ifelse(-μ ==2,1,0)
                iz = iz1 + ifelse(-μ ==3,1,0)
                it = it1 + ifelse(-μ ==4,1,0)

                η = staggered_phase(-μ,ix1,iy1,iz1,it1,NX,NY,NZ,NT)

                ix,iy,iz,it,sign =  apply_periodicity(ix,iy,iz,it,NX,NY,NZ,NT,b.BoundaryCondition)

                outindices[k] = (ix,iy,iz,it,ialpha)

                for k1=1:NC
                    b[k1,ix,iy,iz,it,ialpha] = 0
                    for k2=1:NC
                        b[k1,ix,iy,iz,it,ialpha] += sign*η*conj(u[-μ][k2,k1,ix1,iy1,iz1,it1])*a[k2,ix1,iy1,iz1,it1,ialpha]
                    end
                end

            end

        end

        return outindices

    end

    function fermion_shift!(b::StaggeredFermion,evensite,u::Array{T,1},μ::Int,a::StaggeredFermion) where T <: SU2GaugeFields
        if μ == 0
            substitute!(b,a)
            return
        end

        ibush = ifelse(evensite,0,1)

        NX = a.NX
        NY = a.NY
        NZ = a.NZ
        NT = a.NT
        NC = a.NC

        #NTrange = get_looprange(NT)
        #println(NTrange)
        if μ > 0
            #idel = zeros(Int64,4)
            #idel[μ] = 1

            n6 = size(a.f)[6]
            for ialpha=1:1
                #for it=NTrange
                for it=1:NT
                    it1 = it + ifelse(μ ==4,1,0) #idel[4]
                    for iz=1:NZ
                        iz1 = iz + ifelse(μ ==3,1,0) #idel[3]
                        for iy=1:NY
                            iy1 = iy + ifelse(μ ==2,1,0) #idel[2]
                            xran =1+(1+ibush+iy+iz+it)%2:2:NX
                            for ix in xran
                                ix1 = ix + ifelse(μ ==1,1,0) #idel[1]
                                η = staggered_phase(μ,ix,iy,iz,it,NX,NY,NZ,NT)

                                b[1,ix,iy,iz,it,ialpha] = η*(u[μ][1,1,ix,iy,iz,it]*a[1,ix1,iy1,iz1,it1,ialpha] +
                                                            u[μ][1,2,ix,iy,iz,it]*a[2,ix1,iy1,iz1,it1,ialpha] )

                                b[2,ix,iy,iz,it,ialpha] = η*(u[μ][2,1,ix,iy,iz,it]*a[1,ix1,iy1,iz1,it1,ialpha] +
                                                            u[μ][2,2,ix,iy,iz,it]*a[2,ix1,iy1,iz1,it1,ialpha] )

                                #=
                                for k1=1:3
                                    b[k1,ix,iy,iz,it,ialpha] = 0
                                    for k2=1:3
                                        b[k1,ix,iy,iz,it,ialpha] += u[μ][k1,k2,ix,iy,iz,it]*a[k2,ix1,iy1,iz1,it1,ialpha]
                                    end
                                end
                                =#
                            end
                        end
                    end
                end
            end

        elseif μ < 0
            #idel = zeros(Int64,4)
            #idel[-μ] = 1
            #n6 = size(b.f)[6]
            for ialpha =1:1
                for it=1:NT
                    it1 = it - ifelse(-μ ==4,1,0) #idel[4]
                    for iz=1:NZ
                        iz1 = iz - ifelse(-μ ==3,1,0) #idel[3]
                        for iy=1:NY
                            iy1 = iy - ifelse(-μ ==2,1,0)  #idel[2]
                            xran =1+(1+ibush+iy+iz+it)%2:2:NX
                            for ix in xran
                                ix1 = ix - ifelse(-μ ==1,1,0) #idel[1]
                                η = staggered_phase(-μ,ix1,iy1,iz1,it1,NX,NY,NZ,NT)

                                b[1,ix,iy,iz,it,ialpha] = η*(conj(u[-μ][1,1,ix1,iy1,iz1,it1])*a[1,ix1,iy1,iz1,it1,ialpha] +
                                                            conj(u[-μ][2,1,ix1,iy1,iz1,it1])*a[2,ix1,iy1,iz1,it1,ialpha] )

                                b[2,ix,iy,iz,it,ialpha] = η*(conj(u[-μ][1,2,ix1,iy1,iz1,it1])*a[1,ix1,iy1,iz1,it1,ialpha] +
                                                            conj(u[-μ][2,2,ix1,iy1,iz1,it1])*a[2,ix1,iy1,iz1,it1,ialpha] )

                            end
                        end
                    end
                end
            end
        end

    end

    function fermion_shiftB!(b::StaggeredFermion,u::Array{GaugeFields{SU{NC}},1},μ,a::StaggeredFermion)  where NC
        if μ == 0
            substitute!(b,a)
            return
        end

        NX = a.NX
        NY = a.NY
        NZ = a.NZ
        NT = a.NT
        #NC = a.NC

        if μ > 0
            error("""
            Sorry this case if not yet ready
            mu = $mu
            """)

        elseif μ < 0
            #idel = zeros(Int64,4)
            #idel[-μ] = 1
            #n6 = size(b.f)[6]
            for ialpha =1:1
                for it=1:NT
                    it1 = it + ifelse(-μ == 4,1,0) #idel[4]
                    for iz=1:NZ
                        iz1 = iz +ifelse(-μ == 3,1,0)  #idel[3]
                        for iy=1:NY
                            iy1 = iy + ifelse(-μ == 2,1,0) #idel[2]
                            for ix=1:NX
                                ix1 = ix + ifelse(-μ == 1,1,0)  #idel[1]
                                η = staggered_phase(-μ,ix,iy,iz,it,NX,NY,NZ,NT)

                                for k1=1:NC
                                    b[k1,ix,iy,iz,it,ialpha] = 0
                                    for k2=1:NC
                                        b[k1,ix,iy,iz,it,ialpha] += η*conj(a[k2,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][k1,k2,ix,iy,iz,it])
                                    end
                                end

                            end
                        end
                    end
                end
            end
        end

    end

    function fermion_shiftB!(b::StaggeredFermion,evensite,u::Array{GaugeFields{SU{NC}},1},μ,a::StaggeredFermion)  where NC #T <: Union{SUNGaugeFields,U1GaugeFields}
        if μ == 0
            substitute!(b,a)
            return
        end

        ibush = ifelse(evensite,0,1)

        NX = a.NX
        NY = a.NY
        NZ = a.NZ
        NT = a.NT
        #NC = a.NC

        if μ > 0
            error("""
            Sorry this case if not yet ready
            mu = $mu
            """)

        elseif μ < 0
            #idel = zeros(Int64,4)
            #idel[-μ] = 1
            #n6 = size(b.f)[6]
            for ialpha =1:1
                for it=1:NT
                    it1 = it + ifelse(-μ == 4,1,0) #idel[4]
                    for iz=1:NZ
                        iz1 = iz +ifelse(-μ == 3,1,0)  #idel[3]
                        for iy=1:NY
                            iy1 = iy + ifelse(-μ == 2,1,0) #idel[2]
                            xran =1+(1+ibush+iy+iz+it)%2:2:NX
                            for ix in xran
                            #for ix=1:NX
                                ix1 = ix + ifelse(-μ == 1,1,0)  #idel[1]
                                η = staggered_phase(-μ,ix,iy,iz,it,NX,NY,NZ,NT)

                                for k1=1:NC
                                    b[k1,ix,iy,iz,it,ialpha] = 0
                                    for k2=1:NC
                                        b[k1,ix,iy,iz,it,ialpha] += η*conj(a[k2,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][k1,k2,ix,iy,iz,it])
                                    end
                                end

                            end
                        end
                    end
                end
            end
        end

    end

    function fermion_shiftB!(b::StaggeredFermion,u::Array{SU3GaugeFields,1},μ,a::StaggeredFermion)
        if μ == 0
            substitute!(b,a)
            return
        end

        NX = a.NX
        NY = a.NY
        NZ = a.NZ
        NT = a.NT
        NC = a.NC

        if μ > 0
            error("""
            Sorry this case if not yet ready
            mu = $mu
            """)

        elseif μ < 0
            #idel = zeros(Int64,4)
            #idel[-μ] = 1
            #n6 = size(b.f)[6]
            for ialpha =1:1
                for it=1:NT
                    it1 = it + ifelse(-μ == 4,1,0) #idel[4]
                    for iz=1:NZ
                        iz1 = iz +ifelse(-μ == 3,1,0)  #idel[3]
                        for iy=1:NY
                            iy1 = iy + ifelse(-μ == 2,1,0) #idel[2]
                            for ix=1:NX
                                ix1 = ix + ifelse(-μ == 1,1,0)  #idel[1]
                                η = staggered_phase(-μ,ix,iy,iz,it,NX,NY,NZ,NT)

                                b[1,ix,iy,iz,it,ialpha] = η*(conj(a[1,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][1,1,ix,iy,iz,it]) +
                                                            conj(a[2,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][1,2,ix,iy,iz,it]) +
                                                            conj(a[3,ix1,iy1,iz1,it1,ialpha])* conj(u[-μ][1,3,ix,iy,iz,it]))

                                b[2,ix,iy,iz,it,ialpha] = η*(conj(a[1,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][2,1,ix,iy,iz,it])+
                                                            conj(a[2,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][2,2,ix,iy,iz,it]) +
                                                            conj(a[3,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][2,3,ix,iy,iz,it]))

                                b[3,ix,iy,iz,it,ialpha] = η*(conj(a[1,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][3,1,ix,iy,iz,it]) +
                                                            conj(a[2,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][3,2,ix,iy,iz,it]) +
                                                            conj(a[3,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][3,3,ix,iy,iz,it]))

                                #=
                                for k1=1:NC
                                    b[k1,ix,iy,iz,it,ialpha] = 0
                                    for k2=1:NC
                                        b[k1,ix,iy,iz,it,ialpha] += conj(a[k2,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][k1,k2,ix,iy,iz,it])
                                    end
                                end
                                =#

                            end
                        end
                    end
                end
            end
        end

    end

    function fermion_shiftB!(b::StaggeredFermion,evensite,u::Array{SU3GaugeFields,1},μ,a::StaggeredFermion)
        if μ == 0
            substitute!(b,a)
            return
        end

        ibush = ifelse(evensite,0,1)

        NX = a.NX
        NY = a.NY
        NZ = a.NZ
        NT = a.NT
        NC = a.NC

        if μ > 0
            error("""
            Sorry this case if not yet ready
            mu = $mu
            """)

        elseif μ < 0
            #idel = zeros(Int64,4)
            #idel[-μ] = 1
            #n6 = size(b.f)[6]
            for ialpha =1:1
                for it=1:NT
                    it1 = it + ifelse(-μ == 4,1,0) #idel[4]
                    for iz=1:NZ
                        iz1 = iz +ifelse(-μ == 3,1,0)  #idel[3]
                        for iy=1:NY
                            iy1 = iy + ifelse(-μ == 2,1,0) #idel[2]
                            xran =1+(1+ibush+iy+iz+it)%2:2:NX
                            for ix in xran
                            #for ix=1:NX
                                ix1 = ix + ifelse(-μ == 1,1,0)  #idel[1]
                                η = staggered_phase(-μ,ix,iy,iz,it,NX,NY,NZ,NT)

                                b[1,ix,iy,iz,it,ialpha] = η*(conj(a[1,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][1,1,ix,iy,iz,it]) +
                                                            conj(a[2,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][1,2,ix,iy,iz,it]) +
                                                            conj(a[3,ix1,iy1,iz1,it1,ialpha])* conj(u[-μ][1,3,ix,iy,iz,it]))

                                b[2,ix,iy,iz,it,ialpha] = η*(conj(a[1,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][2,1,ix,iy,iz,it])+
                                                            conj(a[2,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][2,2,ix,iy,iz,it]) +
                                                            conj(a[3,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][2,3,ix,iy,iz,it]))

                                b[3,ix,iy,iz,it,ialpha] = η*(conj(a[1,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][3,1,ix,iy,iz,it]) +
                                                            conj(a[2,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][3,2,ix,iy,iz,it]) +
                                                            conj(a[3,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][3,3,ix,iy,iz,it]))

                                #=
                                for k1=1:NC
                                    b[k1,ix,iy,iz,it,ialpha] = 0
                                    for k2=1:NC
                                        b[k1,ix,iy,iz,it,ialpha] += conj(a[k2,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][k1,k2,ix,iy,iz,it])
                                    end
                                end
                                =#

                            end
                        end
                    end
                end
            end
        end

    end

    function fermion_shiftB!(b::StaggeredFermion,u::Array{SU2GaugeFields,1},μ,a::StaggeredFermion)
        if μ == 0
            substitute!(b,a)
            return
        end

        NX = a.NX
        NY = a.NY
        NZ = a.NZ
        NT = a.NT
        NC = a.NC

        if μ > 0
            error("""
            Sorry this case if not yet ready
            mu = $mu
            """)

        elseif μ < 0
            #idel = zeros(Int64,4)
            #idel[-μ] = 1
            #n6 = size(b.f)[6]
            for ialpha =1:1
                for it=1:NT
                    it1 = it + ifelse(-μ == 4,1,0) #idel[4]
                    for iz=1:NZ
                        iz1 = iz +ifelse(-μ == 3,1,0)  #idel[3]
                        for iy=1:NY
                            iy1 = iy + ifelse(-μ == 2,1,0) #idel[2]
                            for ix=1:NX
                                ix1 = ix + ifelse(-μ == 1,1,0)  #idel[1]
                                η = staggered_phase(-μ,ix,iy,iz,it,NX,NY,NZ,NT)

                                b[1,ix,iy,iz,it,ialpha] = η*(conj(a[1,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][1,1,ix,iy,iz,it]) +
                                                            conj(a[2,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][1,2,ix,iy,iz,it]))

                                b[2,ix,iy,iz,it,ialpha] = η*(conj(a[1,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][2,1,ix,iy,iz,it])+
                                                            conj(a[2,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][2,2,ix,iy,iz,it]) )

                            end
                        end
                    end
                end
            end
        end

    end

    function fermion_shiftB!(b::StaggeredFermion,u::Array{SU2GaugeFields,1},μ,a::StaggeredFermion,indices)
        if μ == 0
            substitute!(b,a)
            return
        end

        NX = a.NX
        NY = a.NY
        NZ = a.NZ
        NT = a.NT
        NC = a.NC

        if μ > 0
            error("""
            Sorry this case if not yet ready
            mu = $mu
            """)

        elseif μ < 0

            ix1 = indices[1]
            iy1 = indices[2]
            iz1 = indices[3]
            it1 = indices[4]
            ialpha = indices[5]

            ix = ix1 - ifelse(-μ ==1,1,0)
            iy = iy1 - ifelse(-μ ==2,1,0)
            iz = iz1 - ifelse(-μ ==3,1,0)
            it = it1 - ifelse(-μ ==4,1,0)

            ix,iy,iz,it,sign =  apply_periodicity(ix,iy,iz,it,NX,NY,NZ,NT,b.BoundaryCondition)

            η = staggered_phase(-μ,ix,iy,iz,it,NX,NY,NZ,NT)

            b[1,ix,iy,iz,it,ialpha] = sign*η*(conj(a[1,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][1,1,ix,iy,iz,it]) +
                                        conj(a[2,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][1,2,ix,iy,iz,it]))

            b[2,ix,iy,iz,it,ialpha] = sign*η*(conj(a[1,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][2,1,ix,iy,iz,it])+
                                        conj(a[2,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][2,2,ix,iy,iz,it]) )

        end

        return (ix,iy,iz,it,ialpha)

    end

    function fermion_shiftB!(b::StaggeredFermion,evensite,u::Array{SU2GaugeFields,1},μ,a::StaggeredFermion)
        if μ == 0
            substitute!(b,a)
            return
        end

        ibush = ifelse(evensite,0,1)

        NX = a.NX
        NY = a.NY
        NZ = a.NZ
        NT = a.NT
        NC = a.NC

        if μ > 0
            error("""
            Sorry this case if not yet ready
            mu = $mu
            """)

        elseif μ < 0
            #idel = zeros(Int64,4)
            #idel[-μ] = 1
            #n6 = size(b.f)[6]
            for ialpha =1:1
                for it=1:NT
                    it1 = it + ifelse(-μ == 4,1,0) #idel[4]
                    for iz=1:NZ
                        iz1 = iz +ifelse(-μ == 3,1,0)  #idel[3]
                        for iy=1:NY
                            iy1 = iy + ifelse(-μ == 2,1,0) #idel[2]
                            xran =1+(1+ibush+iy+iz+it)%2:2:NX
                            for ix in xran
                                ix1 = ix + ifelse(-μ == 1,1,0)  #idel[1]
                                η = staggered_phase(-μ,ix,iy,iz,it,NX,NY,NZ,NT)

                                b[1,ix,iy,iz,it,ialpha] = η*(conj(a[1,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][1,1,ix,iy,iz,it]) +
                                                            conj(a[2,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][1,2,ix,iy,iz,it]))

                                b[2,ix,iy,iz,it,ialpha] = η*(conj(a[1,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][2,1,ix,iy,iz,it])+
                                                            conj(a[2,ix1,iy1,iz1,it1,ialpha])*conj(u[-μ][2,2,ix,iy,iz,it]) )

                            end
                        end
                    end
                end
            end
        end

    end

end