# Measure the kinetic (one-body) and the potential (two-body) energy
function measure_energy(
    lattice::Lattice, UUp_Array::Array{T1,3}, UDn_Array::Array{T2,3},
    n̄_up::Array{T3,1}, n̄_dn::Array{T4,1}
) where {T1<:FloatType, T2<:FloatType, T3<:FloatType, T4<:FloatType}

    EOneBody, ETwoBody = 0, 0

    for i in eachindex(view(lattice.KMatrix, 1 : lattice.Ns, 1 : lattice.Ns))
        if lattice.KMatrix[i] != 0
            Dij_up = sum(UUp_Array[:, i[1], i[2]] .* n̄_up)
            Dij_dn = sum(UDn_Array[:, i[1], i[2]] .* n̄_dn)
            EOneBody += -lattice.t * (Dij_up + Dij_dn)
        end
    end

    for i = 1 : lattice.Ns
        Umatrix = UUp_Array[:, i, i] * UDn_Array[:, i, i]'
        nmatrix = n̄_up * n̄_dn'
        Diiii = sum(Umatrix .* nmatrix)
        ETwoBody += lattice.U * Diiii
    end

    real(EOneBody), real(ETwoBody)
end

function mc_integration(lattice::Lattice, qmc::QMC)

    Ek_Array = Vector{Float64}()    # kinetic energy
    Ep_Array = Vector{Float64}()    # potential energy
    E_Array = Vector{Float64}()     # total energy
    sign_Array = Vector{Float64}()  # sign
    wrapping_count = 0

    # initialize a random field configuration
    σ_field = 2 * (rand(Float64, (lattice.Ns, qmc.L)) .< 0.5) .- 1
    # initial propagation
    Q0_Up, D0_Up, T0_Up, Q0_Dn, D0_Dn, T0_Dn = 
                                    matrix_product_stablized(σ_field, lattice, qmc)

    eigvals_up = eigvals(D0_Up * T0_Up * Q0_Up, sortby = abs)
    eigvals_dn = eigvals(D0_Dn * T0_Dn * Q0_Dn, sortby = abs)
    # recursion
    ZUp = main_recursion(lattice.Ns, lattice.nup, eigvals_up, false)
    ZDn = main_recursion(lattice.Ns, lattice.ndn, eigvals_dn, false)
    # projection
    #ZUp = particle_projection(lattice.Ns, lattice.nup, eigvals_up, false)
    #ZDn = particle_projection(lattice.Ns, lattice.ndn, eigvals_dn, false)

    for n = 1 : qmc.sample_size
        
        # sweep over the field
        for l = 1 : qmc.L
            # recalculate the decomposition over a certain period
            # to ensure numerical stability
            if wrapping_count == 15
                σ_field_shifted = circshift(σ_field, (0, -l + 1))
                Q0_Up, D0_Up, T0_Up, Q0_Dn, D0_Dn, T0_Dn = 
                                    matrix_product_stablized(σ_field_shifted, lattice, qmc)
                wrapping_count = 0
            end
            # store the decomposition temporarily
            Q0_Up_temp, D0_Up_temp, T0_Up_temp = Q0_Up, D0_Up, T0_Up
            Q0_Dn_temp, D0_Dn_temp, T0_Dn_temp = Q0_Dn, D0_Dn, T0_Dn
            # (use the first time slice as an example)
            # the following lines convert B_L * ... * B_1 = Q_0 * D_0 * T_0 to
            # B_1 * B_L * ... * B_2 = B_1 * [Q_0 * D_0 * T_0 * B_1^-1] = B_1 * Q_l * D_l * T_l
            # which preserves the eigenvalues and the trace.
            σ = σ_field[:, l]
            Bl_Up, Bl_Dn = singlestep_matrix(σ, lattice, qmc)

            QRCP_update!(Q0_Up, D0_Up, T0_Up, inv(Bl_Up), 'R')
            QRCP_update!(Q0_Dn, D0_Dn, T0_Dn, inv(Bl_Dn), 'R')

            # loop over sites
            for i = 1 : lattice.Ns

                # propose a new configuration
                σ[i] *= -1
                Bl_Up, Bl_Dn = singlestep_matrix(σ, lattice, qmc)

                Ql_Up, Dl_Up, Tl_Up = QRCP_update!(copy(Q0_Up), copy(D0_Up), copy(T0_Up),
                                                    Bl_Up, 'L')
                Ql_Dn, Dl_Dn, Tl_Dn = QRCP_update!(copy(Q0_Dn), copy(D0_Dn), copy(T0_Dn),
                                                    Bl_Dn, 'L')
                
                eigvals_up = eigvals(Dl_Up * Tl_Up * Ql_Up, sortby = abs)
                eigvals_dn = eigvals(Dl_Dn * Tl_Dn * Ql_Dn, sortby = abs)
                # recursion
                ZUpNew = main_recursion(lattice.Ns, lattice.nup, eigvals_up, false)
                ZDnNew = main_recursion(lattice.Ns, lattice.ndn, eigvals_dn, false)
                # projection
                #ZUpNew = particle_projection(lattice.Ns, lattice.nup, eigvals_up, false)
                #ZDnNew = particle_projection(lattice.Ns, lattice.ndn, eigvals_dn, false)
                # caculate the acceptence ratio
                R = abs((ZUpNew / ZUp) * (ZDnNew / ZDn))
                accept_ratio = R / (1 + R)

                if rand() < accept_ratio

                    # accept, update the field
                    σ_field[i, l] *= -1
                    ZUp, ZDn = ZUpNew, ZDnNew
                    # store the decomposition temporarily
                    Q0_Up_temp, D0_Up_temp, T0_Up_temp = Ql_Up, Dl_Up, Tl_Up
                    Q0_Dn_temp, D0_Dn_temp, T0_Dn_temp = Ql_Dn, Dl_Dn, Tl_Dn

                else
                    # roll back the change
                    σ[i] *= -1
                end

                # update Q_0, D_0, T_0 at the last site
                if i == lattice.Ns
                    Q0_Up, D0_Up, T0_Up = Q0_Up_temp, D0_Up_temp, T0_Up_temp
                    Q0_Dn, D0_Dn, T0_Dn = Q0_Dn_temp, D0_Dn_temp, T0_Dn_temp
                end

            end

            wrapping_count += 1

        end

        # measure observables
        # propagate using normal order
        MatProdUp_Q, MatProdUp_D, MatProdUp_T, MatProdDn_Q, MatProdDn_D, MatProdDn_T =
                                        matrix_product_stablized(σ_field, lattice, qmc)
        eigen_up = eigen(MatProdUp_D * MatProdUp_T * MatProdUp_Q, sortby = abs)
        eigen_dn = eigen(MatProdDn_D * MatProdDn_T * MatProdDn_Q, sortby = abs)
        # transformation matrix (A = P D P^-1)
        MatProdUp_P = MatProdUp_Q * eigen_up.vectors
        MatProdDn_P = MatProdDn_Q * eigen_dn.vectors
        # overlap coefficient
        UUp_Array, UDn_Array = overlap_coefficient(lattice.Ns, MatProdUp_P, MatProdDn_P)
        # recursion
        ZUp, n̄up = main_recursion(lattice.Ns, lattice.nup, eigen_up.values, true)
        ZDn, n̄dn = main_recursion(lattice.Ns, lattice.ndn, eigen_dn.values, true)
        # projection
        #ZUp, n̄up = particle_projection(lattice.Ns, lattice.nup, eigen_up.values, true)
        #ZDn, n̄dn = particle_projection(lattice.Ns, lattice.ndn, eigen_dn.values, true)
        # sign
        sgn = sign(ZUp) * sign(ZDn)
        # measure energy
        Ek, Ep = measure_energy(lattice, UUp_Array, UDn_Array, n̄up, n̄dn)
        # store data
        push!(sign_Array, sgn)
        push!(Ek_Array, Ek * sgn)
        push!(Ep_Array, Ep * sgn)
        push!(E_Array, (Ek + Ep) * sgn)

    end

    # block average
    EMean = mean(E_Array)
    EVar = varm(E_Array, EMean) / qmc.sample_size
    EkMean = mean(Ek_Array)
    EkVar = varm(Ek_Array, EkMean) / qmc.sample_size
    EpMean = mean(Ep_Array)
    EpVar = varm(Ep_Array, EpMean) / qmc.sample_size
    avgsgn = mean(sign_Array)

    EMean, EVar, EkMean, EkVar, EpMean, EpVar, avgsgn
    
end
