text:
module PeriodizeSC

using Cubature: hcubature
using Cuba
using JSON
using NPZ

_Nc = 4
_NcN = 8 #Sites for nambu = _Nc*2

# Global Variables
II = eye(Complex{Float64}, _Nc)
ZEROS = zeros(Complex{Float64}, _Nc, _Nc)

"""
Type ModelVector holds all the relevant quantities needed for computations. The built-in member function ModelVector() creates an
instance of the type ModelVector.

#Arguments of ModelVector():

- t_::Float64: Nearest-neighbor hopping amplitude (Float64)
- tp_::Float64: Second nearest-neighbor hopping amplitude (Float64)
- tpp_::Float64: Third nearest-neighbor hopping amplitude (Float64)
- mu_::Float64: Chemical potential (Float64)
- wvec_::Float64: Matsubara frequency grid (Array{Complex{Float64}, 1})
- sEvec_c::Array{Complex{Float64}, 3}: Vector of self-energies (Array{Complex{Float64}, 3})
- cumulants_::Array{Complex{Float64}, 3}: Vector of cumulants of the Green's function (Array{Complex{Float64}, 3})

#Returns:

- ModelVector instance.

"""
type ModelVector
   t_::Float64 ; tp_::Float64 ; tpp_::Float64;  mu_::Float64
   wvec_::Array{Complex{Float64}, 1} ; sEvec_c_::Array{Complex{Float64}, 3}
   cumulants_::Array{Complex{Float64}, 3}

   function ModelVector(t::Float64, tp::Float64, tpp::Float64, mu::Float64,
                        wvec::Array{Complex{Float64}, 1}, sEvec_c::Array{Complex{Float64}, 3})

        cumulants = build_cumulants(wvec, mu, sEvec_c)
        return new(t, tp, tpp, mu, wvec, sEvec_c, cumulants)
    end

end

"""
Function member of the type ModelVector returning the cumulant of the Green's function.
Useful to compute the superfluid stiffness by periodizing the cumulant.

#Arguments:

- wvec_::Array{Complex{Float64}, 1}: Matsubara frequency grid (Array{Complex{Float64}, 1})
- mu_::Float64: Chemical potential (Float64)
- sEvec_c_::Array{Complex{Float64}, 3}: Vector of self-energies (Array{Complex{Float64}, 3})

#Returns:

- cumulants: Cumulants of the Green's function (Complex{Float64} matrix of size 2*Nc X 2*Nc)

"""
function build_cumulants(wvec::Array{Complex{Float64}, 1}, mu::Float64, sEvec_c::Array{Complex{Float64}, 3})

      cumulants = zeros(Complex{Float64}, size(sEvec_c))

      for (ii, ww) in enumerate(wvec)
          tmp = zeros(Complex{Float64}, (_NcN, _NcN))
          tmp[1, 1] = tmp[2, 2] = tmp[3, 3] = tmp[4, 4] = (ww + mu)
          tmp[5, 5] = tmp[6, 6] = tmp[7, 7] = tmp[8, 8] = -conj((ww + mu))
          tmp -= sEvec_c[ii, :, :]
          cumulants[ii, :, :] = inv(tmp)
      end

      return cumulants
end

"""
Type Model holds all the relevant quantities needed for computations for a given Matsubara frequency. The built-in member function Model() instantiates
an object that a priori has a type ModelVector.

#Arguments of Model():

- modelvec: Object of type ModelVector
- ii::Integer: Matsubara frequency index number (Int64)

#Returns:

- Model instance.

"""
type Model
    t_::Float64  ; tp_::Float64 ; tpp_::Float64 ; mu_::Float64
    w_::Complex{Float64} ; sE_::Array{Complex{Float64}, 2}
    cumulant_::Array{Complex{Float64}, 2}

    function Model(modelvec::ModelVector, ii::Integer)
        (t, tp, tpp, mu, w, sE_c, cumulant) = (modelvec.t_, modelvec.tp_, modelvec.tpp_, modelvec.mu_, modelvec.wvec_[ii], modelvec.sEvec_c_[ii, :, :], modelvec.cumulants_[ii, :, :])
        return new(t, tp, tpp, mu, w, sE_c, cumulant)
    end

end

"""
Function to construct the nearest-neighbor dispersion relation on CuO2 square lattice.

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction

#Returns:

- epsilonk: Value of the dispersion relation on the 2d reciprocal lattice

"""
function epsilonk(model::Model, kx::Float64, ky::Float64)
    epsilonk = -2.0*model.t_*(cos(kx) + cos(ky))
    return epsilonk
end

"""
Function to construct the nearest-neighbor dispersion relation derivative with respect to kx on CuO2 square lattice.

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction

#Returns:

- dxepsilonk: Value of the dispersion relation derivative along kx on the 2d reciprocal lattice

"""
function Dxepsilonk(model::Model, kx::Float64, ky::Float64)
    dxepsilonk = 2.0*model.t_*sin(kx)
    return dxepsilonk
end

"""
Function to construct the nearest-neighbor dispersion relation derivative with respect to ky on CuO2 square lattice.

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction

#Returns:

- dyepsilonk: Value of the dispersion relation derivative along ky on the 2d reciprocal lattice

"""
function Dyepsilonk(model::Model, kx::Float64, ky::Float64)
    dyepsilonk = 2.0*model.t_*sin(ky)
    return dyepsilonk
end

"""
Function to construct the second nearest-neighbor dispersion relation on CuO2 square lattice.

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction

#Returns:

- zetak: Value of the dispersion relation on the 2d reciprocal lattice

"""
function zetak(model::Model, kx::Float64, ky::Float64)
    zetak = -2.0*model.tp_*(cos(kx+ky) + cos(kx-ky))
    return zetak
end

"""
Function to construct the second nearest-neighbor dispersion relation derivative along kx on CuO2 square lattice.

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction

#Returns:

- dxzetak: Value of the dispersion relation derivative along kx on the 2d reciprocal lattice

"""
function Dxzetak(model::Model, kx::Float64, ky::Float64)
    dxzetak = 2.0*model.tp_*(sin(kx+ky)+sin(kx-ky))
    return dxzetak
end

"""
Function to construct the second nearest-neighbor dispersion relation derivative along ky on CuO2 square lattice.

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction

#Returns:

- dyzetak: Value of the dispersion relation derivative along ky on the 2d reciprocal lattice

"""
function Dyzetak(model::Model, kx::Float64, ky::Float64)
    dyzetak = 2.0*model.tp_*(sin(kx+ky)-sin(kx-ky))
    return dyzetak
end

"""
Function to construct the third nearest-neighbor dispersion relation on CuO2 square lattice.

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction

#Returns:

- omegak: Value of the dispersion relation on the 2d reciprocal lattice

"""
function omegak(model::Model, kx::Float64, ky::Float64)
    omegak = -2.0*model.tpp_*(cos(2.0*kx)+cos(2.0*ky))
    return omegak
end

"""
Function to construct the third nearest-neighbor dispersion relation derivative along kx on CuO2 square lattice

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction

#Returns:

- dxomegak: Value of the dispersion relation derivative along kx on the 2d reciprocal lattice

"""
function Dxomegak(model::Model, kx::Float64, ky::Float64)
    dxomegak = 4.0*model.tpp_*sin(2.0*kx)
    return dxomegak
end

"""
Function to construct the third-nearest-neighbor dispersion relation derivative along ky on CuO2 square lattice

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction

#Returns:

- dyomegak: Value of the dispersion relation derivative along ky on the 2d reciprocal lattice

"""
function Dyomegak(model::Model, kx::Float64, ky::Float64)
    dyomegak = 4.0*model.tpp_*sin(2.0*ky)
    return dyomegak
end

"""
Function to construct the full dispersion relation on CuO2 square lattice

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction

#Returns:

- eps_0: Value of the full dispersion relation on the 2d reciprocal lattice

"""
function eps_0(model::Model, kx::Float64, ky::Float64)
    return(epsilonk(model,kx,ky)+zetak(model,kx,ky)+omegak(model,kx,ky))
end

############################################################################## Important Note ###########################################################################
## The following functions construct the bare current vertices used when computing the in-plane superfluid stiffness. The situation in which each and every one of the ##
## functions is used is specified.                                                                                ##
#########################################################################################################################################################################

######## Cumulant and Green's function periodization schemes SC only
function DxDxEpsilonbark(model::Model, kx::Float64, ky::Float64)
    return (Dxepsilonk(model,kx,ky)+Dxzetak(model,kx,ky)+Dxomegak(model,kx,ky))*(Dxepsilonk(model,kx,ky)+Dxzetak(model,kx,ky)+Dxomegak(model,kx,ky))
end

function DyDyEpsilonbark(model::Model, kx::Float64, ky::Float64)
    return (Dyepsilonk(model,kx,ky)+Dyzetak(model,kx,ky)+Dyomegak(model,kx,ky))*(Dyepsilonk(model,kx,ky)+Dyzetak(model,kx,ky)+Dyomegak(model,kx,ky))
end

####### Cumulant and Green's function periodization schemes SC+AFM mixed state
function DxDxZetakbar(model::Model, kx::Float64, ky::Float64)
    return (Dxzetak(model,kx,ky)+Dxomegak(model,kx,ky))*(Dxzetak(model,kx,ky)+Dxomegak(model,kx,ky))
end

function DyDyZetakbar(model::Model, kx::Float64, ky::Float64)
    return (Dyzetak(model,kx,ky)+Dyomegak(model,kx,ky))*(Dyzetak(model,kx,ky)+Dyomegak(model,kx,ky))
end

function DxDxZetakbarepsilonk(model::Model, kx::Float64, ky::Float64)
    return (Dxzetak(model,kx,ky)+Dxomegak(model,kx,ky))*Dxepsilonk(model,kx,ky)
end

function DyDyZetakbarepsilonk(model::Model, kx::Float64, ky::Float64)
    return (Dyzetak(model,kx,ky)+Dyomegak(model,kx,ky))*Dyepsilonk(model,kx,ky)
end

######## Trace scheme for SC+AFM mixed state
function DxDxzetak(model::Model, kx::Float64, ky::Float64)
    return Dxzetak(model,kx,ky)*Dxzetak(model,kx,ky)
end

function DyDyzetak(model::Model, kx::Float64, ky::Float64)
    return Dyzetak(model,kx,ky)*Dyzetak(model,kx,ky)
end

function DxDxomegak(model::Model, kx::Float64, ky::Float64)
    return Dxomegak(model,kx,ky)*Dxomegak(model,kx,ky)
end

function DyDyomegak(model::Model, kx::Float64, ky::Float64)
    return Dyomegak(model,kx,ky)*Dyomegak(model,kx,ky)
end

function DxDxepsilonk(model::Model, kx::Float64, ky::Float64)
    return Dxepsilonk(model,kx,ky)*Dxepsilonk(model,kx,ky)
end

function DyDyepsilonk(model::Model, kx::Float64, ky::Float64)
    return Dyepsilonk(model,kx,ky)*Dyepsilonk(model,kx,ky)
end

function DxDxzetakomegak(model::Model, kx::Float64, ky::Float64)
    return Dxzetak(model,kx,ky)*Dxomegak(model,kx,ky)
end

function DyDyzetakomegak(model::Model, kx::Float64, ky::Float64)
    return Dyzetak(model,kx,ky)*Dyomegak(model,kx,ky)
end

function DxDxzetakepsilonk(model::Model, kx::Float64, ky::Float64)
    return Dxzetak(model,kx,ky)*Dxepsilonk(model,kx,ky)
end

function DyDyzetakepsilonk(model::Model, kx::Float64, ky::Float64)
    return Dyzetak(model,kx,ky)*Dyepsilonk(model,kx,ky)
end

function DxDxomegakepsilonk(model::Model, kx::Float64, ky::Float64)
    return Dxomegak(model,kx,ky)*Dxepsilonk(model,kx,ky)
end

function DyDyomegakepsilonk(model::Model, kx::Float64, ky::Float64)
    return Dyomegak(model,kx,ky)*Dyepsilonk(model,kx,ky)
end

###################################################################### End of Important Note ################################################################################
#############################################################################################################################################################################

"""
Function defining the hopping between two layers (CuO2 planes) within a unit cell.

#Arguments:

- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction

#Returns:

- tperp_squared: Value of the bilayer hopping in terms of the in-plane wavevectors
(Represents the current vertices in the current-current correlation function)

"""
function tperp(kx::Float64,ky::Float64)
    coskx = cos(kx)
    cosky = cos(ky)
    tperp = -(coskx-cosky)*(coskx-cosky)
    tperp_squared = tperp*tperp
    return 1/2*tperp_squared  ### In reality, one must add a factor of 1/2 to account of the integral along z-axis.
end

"""
Function adding the proper phases proportional to the reduced supercluster BZ to the intercluster hoppings.

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction

#Returns:

- 2N_c x 2N_c complex-valued matrix t(k^{tilde})
"""
function tktilde(model::Model, kx::Float64, ky::Float64)
    t = model.t_ ; tp = model.tp_; tpp = model.tpp_
    k = [kx, ky]
    r_sites = [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]
    K_sites = pi*deepcopy(r_sites)
    t_array = zeros(Complex{Float64}, (_Nc, _Nc))

    for i in 1:_Nc
        for j in 1:_Nc
            for K in K_sites
                t_array[i, j] += 1.0/_Nc * exp(1.0im*dot(K + k, r_sites[i] - r_sites[j])) * eps_0(model, (K + k)...)
            end
        end
    end

    return (vcat(hcat(t_array, ZEROS), hcat(ZEROS, -t_array)))
end

"""
Function building the full cluster Green's function.

#Arguments:

- model: model instance to have access to the attributes such as hopping elements, self-energy and chemical potential
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction

#Returns:

- gf_ktilde: 2N_c x 2N_c complex-valued matrix G_c(i omega_n,k^{tilde})
"""
function build_gf_ktilde(model::Model, kx::Float64, ky::Float64)
    zz = (model.w_ + model.mu_)*II
    zz_c = -conj(-conj(model.w_) + model.mu_)*II
    tmp_zz = (vcat(hcat(zz, ZEROS), hcat(ZEROS, -conj(zz)))) #### -conj(zz)
    gf_ktilde = inv(tmp_zz - tktilde(model, kx, ky) - model.sE_)
    return gf_ktilde
end

"""
Function pre-computing the full cluster Green's function on the BZ. Useful to speed up in-plane computations.

#Arguments:

- model: model instance to have access to the attributes such as hopping elements, self-energy and chemical potential
- tk_el::Array{Complex{Float64},2}: pre-computed k-tilde-dependent hopping term (t(k^{tilde}))
#Returns:

- gf_ktilde: 2N_c x 2N_c complex-valued matrix G_c(i omega_n,k^{tilde})
"""
function build_gf_ktilde_prebuild(model::Model, tk_el::Array{Complex{Float64},2})
    zz = (model.w_ + model.mu_)*II
    tmp_zz = (vcat(hcat(zz, ZEROS), hcat(ZEROS, -conj(zz))))  #Same thing as zz_c = -conj(-conj(model.w_) + model.mu_)*II
    gf_ktilde = inv(tmp_zz - tk_el - model.sE_)
    return gf_ktilde
end

"""
Function used to periodize the Green's function with SC only.

#Arguments:

- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction
- args::Array{Complex{Float64},2}: Lattice Green's function in terms of k_tilde before being periodized (Array{Complex{Float64}, 2})

#Returns:

- nambu_periodized: Periodized lattice Green's function (Complex{Float64}-valued matrix 0.5*Nc X 0.5*Nc)

"""
function periodize_nocoex(kx::Float64, ky::Float64, arg::Array{Complex{Float64},2})
    ex = exp(1.0im*kx)
    ey = exp(1.0im*ky)
    emx = conj(ex)
    emy = conj(ey)
    vk1 = [1.0,-1.0*ex,-1.0*ey,ex*ey]                      #equivalent transformation matrix for Green per. would be:
    vk2 = [-1.0*emx,1.0,emx*ey,-1.0*ey]                    #[1.0,-1.0*ex,-1.0*ey,ex*ey]
    vk3 = [-1.0*emy,emy*ex,1.0,-1.0*ex]                    #[-1.0*emx,1.0,emx*ey,-1.0*ey]
    vk4 = [emx*emy,-1.0*emy,-1.0*emx,1.0]                  #[-1.0*emy,emy*ex,1.0,-1.0*ex]
    nambu_periodized = zeros(Complex{Float64},(2,2))       #[emx*emy,-1.0*emy,-1.0*emx,1.0]

    gup = arg[1:4,1:4]
    ff = arg[1:4,5:end]
    ffdag = arg[5:end,1:4]
    gdown = arg[5:end,5:end]

    llperiodized = [gup,ff,ffdag,gdown]
    vk = [vk1,vk2,vk3,vk4]

    list_total = []
    for elem in llperiodized
        list_i = []
        summ=0.0
        for i in 1:size(elem)[1]
            for j in 1:size(elem)[2]
                push!(list_i, elem[i,j]*vk[i][j])
            end
        end
        summ=sum(list_i)
        #println("sum = ", summ)
        push!(list_total,summ)
    end

    nambu_periodized[1,1] = list_total[1]
    nambu_periodized[1,2] = list_total[2]
    nambu_periodized[2,1] = list_total[3]
    nambu_periodized[2,2] = list_total[4]

    return 0.25*nambu_periodized
end

"""
Function used to periodize the cumulant with SC only.

#Arguments:

- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction
- K::Array{Float64,1}: Array of superlattice reciprocal wavevectors (noted K)
- args::Array{Complex{Float64},2}: Lattice Green's function in terms of k_tilde before being periodized (Array{Complex{Float64}, 2})

#Returns:

- nambu_cum: Periodized lattice Green's function (Complex{Float64}-valued matrix 0.5*Nc X 0.5*Nc)

"""
function periodize_nocoex_cum(kx::Float64, ky::Float64, K::Array{Float64,1}, arg::Array{Complex{Float64},2})
    nambu_cum = zeros(Complex{Float64},(2,2))
    R = [[0,0],[1,0],[0,1],[1,1]]
    k = [kx,ky]
    gup = arg[1:4,1:4]
    ff = arg[1:4,5:end]
    ffdag = arg[5:end,1:4]
    gdown = arg[5:end,5:end]

    llperiodized = [gup,ff,ffdag,gdown]
    nambu_sum = Array{Complex{Float64},1}()
    for elem in llperiodized
        tot = 0.0
        for (i,R1) in enumerate(R)
            for (j,R2) in enumerate(R)
                tot += exp(-1.0im*dot(k + K + [pi,pi],R1 - R2))*elem[i,j]
            end
        end
        push!(nambu_sum,tot)
    end
    nambu_cum[1,1] = nambu_sum[1]; nambu_cum[1,2] = nambu_sum[2]
    nambu_cum[2,1] = nambu_sum[3]; nambu_cum[2,2] = nambu_sum[4]

    return nambu_cum
end

"""
Function associated to the function periodizing the cumulant with SC only. Adds the dispersion relations after having built the cumulant.

#Arguments:

- model: model instance to have access to the attributes such as hopping elements, self-energy and chemical potential
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction
- K::Array{Float64,1}: Array of superlattice reciprocal wavevectors (noted K)

#Returns:

- Cumulant of the Green's function (Complex{Float64}-valued matrix Nc X Nc)

"""
function periodize_nocoex_cum_suite(model::Model, kx::Float64, ky::Float64, K::Array{Float64,1})
    nambu_cum_inv = zeros(Complex{Float64},(2,2))

    nambu_cum_inv = inv(periodize_nocoex_cum(kx,ky,K,model.cumulant_))

    k = [kx,ky]

    nambu_cum_inv[1,1] -= eps_0(model,(k + K + [pi,pi])...); nambu_cum_inv[2,2] += eps_0(model,(k + K + [pi,pi])...)

    return inv(nambu_cum_inv)
end

"""
Function associated to the function periodizing the cumulant with SC only. Last inversion step corresping to the
full Green's function.

#Arguments:

- model: model instance to have access to the attributes such as hopping elements, self-energy and chemical potential
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction

#Returns:

- nambu_periodized: Periodized lattice Green's function (Complex{Float64}-valued matrix Nc X Nc)
"""
function periodize_nocoex_cum_finale(model::Model, kx::Float64, ky::Float64)
    nambu_periodized = zeros(Complex{Float64},(2,2))
    Ks = [[0.,0.],[pi,0.],[0.,pi],[pi,pi]]

    for K in Ks
        #println(K)
        nambu_periodized += periodize_nocoex_cum_suite(model,kx,ky,K)
    end

    return 0.25*nambu_periodized
end

"""
Function to periodize the Green's function with coexistence between AF and SC.

#Arguments:

- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction
- args::Array{Complex{Float64},2}: Lattice Green's function in terms of k_tilde before periodized (Array{Complex{Float64}, 2})

#Returns:

- nambu_periodized: Periodized lattice Green's function (Complex{Float64}-valued matrix Nc X Nc)

"""
function periodize_AFM_orb(arg::Array{Complex{Float64},2}, kx::Float64, ky::Float64) ###Appropriate periodization method for AF-SC systems
    R_A = [[0.0,0.0],[1.0,1.0]]
    R_B = [[1.0,0.0],[0.0,1.0]]
    K_x = [[0.0,0.0],[pi,0.0]]; K_y = [[0.0,0.0],[0.0,pi]]
    k = [kx,ky]
    nambu_periodized = zeros(Complex{Float64}, _Nc, _Nc)

    gup = arg[1:4, 1:4]
    ff = arg[1:4, 5:end]
    ffdag = arg[5:end, 1:4]
    gdown = arg[5:end, 5:end]

    llgreen = [gup, ff, ffdag, gdown]
    blocks = zeros(Complex{Float64},2,2,4)
    for (ii,elem) in enumerate(llgreen)
        gAA = [[elem[1,1],elem[1,4]],[elem[4,1],elem[4,4]]]
        gAB = [[elem[1,2],elem[1,3]],[elem[4,2],elem[4,3]]]
        gBA = [[elem[2,1],elem[2,4]],[elem[3,1],elem[3,4]]]
        gBB = [[elem[2,2],elem[2,3]],[elem[3,2],elem[3,3]]]
        summAA = summAB = summBA = summBB = 0.0
        for K in K_y    ####K_x and K_y should be equivalent
            for RA in 1:size(R_A)[1]
                for RAprime in 1:size(R_A)[1]
                    summAA += exp(1.0im*dot(K + k, R_A[RAprime] - R_A[RA]))*gAA[RA][RAprime]
                end
            end
        end
        blocks[1,1,ii] = summAA
        for K in K_y
            for RA in 1:size(R_A)[1]
                for RBprime in 1:size(R_B)[1]
                    summAB += exp(1.0im*dot(K + k, R_B[RBprime] - R_A[RA]))*gAB[RA][RBprime]
                end
            end
        end
        blocks[1,2,ii] = summAB
        for K in K_y
            for RB in 1:size(R_B)[1]
                for RAprime in 1:size(R_A)[1]
                    summBA += exp(1.0im*dot(K + k, R_A[RAprime] - R_B[RB]))*gBA[RB][RAprime]
                end
            end
        end
        blocks[2,1,ii] = summBA
        for K in K_y
            for RB in 1:size(R_B)[1]
                for RBprime in 1:size(R_B)[1]
                    summBB += exp(1.0im*dot(K + k, R_B[RBprime] - R_B[RB]))*gBB[RB][RBprime]
                end
            end
        end
        blocks[2,2,ii] = summBB
    end
    nambu_periodized[1:2,1:2] = blocks[:,:,1]
    nambu_periodized[1:2,3:end] = blocks[:,:,2]
    nambu_periodized[3:end,1:2] = blocks[:,:,3]
    nambu_periodized[3:end,3:end] = blocks[:,:,4]

    return(0.5*nambu_periodized)
end

"""
Function to periodize the cumulant in the case of coexistence between AF and SC.

#Arguments:

- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction
- args::Array{Complex{Float64},2}: Lattice Green's function in terms of k_tilde before periodized (Array{Complex{Float64}, 2})
- K::Array{Float64,1}: Array of superlattice reciprocal wavevectors (noted K)

#Returns:

- nambu_cum: Periodized lattice Green's function (Complex{Float64}-valued matrix of size Nc x Nc)

"""
function periodize_AFM_orb_cum(arg::Array{Complex{Float64},2}, kx::Float64, ky::Float64, K::Array{Float64,1}) ###Appropriate periodization method for AFM-SC systems
    nambu_cum = zeros(Complex{Float64},4,4)
    R_A = [[0.0,0.0],[1.0,1.0]]
    R_B = [[1.0,0.0],[0.0,1.0]]
    k = [kx,ky]

    gup = arg[1:4, 1:4]
    ff = arg[1:4, 5:end]
    ffdag = arg[5:end, 1:4]
    gdown = arg[5:end, 5:end]

    llgreen = [gup, ff, ffdag, gdown]
    sum_blocks = zeros(Complex{Float64},2,2,4)
    for (ii,elem) in enumerate(llgreen)
        gAA = [[elem[1,1],elem[1,4]],[elem[4,1],elem[4,4]]]
        gAB = [[elem[1,2],elem[1,3]],[elem[4,2],elem[4,3]]]
        gBA = [[elem[2,1],elem[2,4]],[elem[3,1],elem[3,4]]]
        gBB = [[elem[2,2],elem[2,3]],[elem[3,2],elem[3,3]]]
        summAA = summAB = summBA = summBB = 0.0
        for RA in 1:size(R_A)[1]
            for RAprime in 1:size(R_A)[1]
                summAA += exp(1.0im*dot(K + k, R_A[RAprime] - R_A[RA]))*gAA[RA][RAprime]
            end
        end
        sum_blocks[1,1,ii] = summAA
        for RA in 1:size(R_A)[1]
            for RBprime in 1:size(R_B)[1]
                summAB += exp(1.0im*dot(K + k, R_B[RBprime] - R_A[RA]))*gAB[RA][RBprime]
            end
        end
        sum_blocks[1,2,ii] = summAB
        for RB in 1:size(R_B)[1]
            for RAprime in 1:size(R_A)[1]
                summBA += exp(1.0im*dot(K + k, R_A[RAprime] - R_B[RB]))*gBA[RB][RAprime]
            end
        end
        sum_blocks[2,1,ii] = summBA
        for RB in 1:size(R_B)[1]
            for RBprime in 1:size(R_B)[1]
                summBB += exp(1.0im*dot(K + k, R_B[RBprime] - R_B[RB]))*gBB[RB][RBprime]
            end
        end
        sum_blocks[2,2,ii] = summBB
    end
    nambu_cum[1:2,1:2] = sum_blocks[:,:,1]
    nambu_cum[1:2,3:end] = sum_blocks[:,:,2]
    nambu_cum[3:end,1:2] = sum_blocks[:,:,3]
    nambu_cum[3:end,3:end] = sum_blocks[:,:,4]

    return nambu_cum
end

"""
Function associated to the function periodizing the cumulant in the case of coexistence between AFM and SC. Adds the
dispersion relations after having built the cumulant.

#Arguments:

- model: model instance to have access to the attributes such as hopping elements, self-energy and chemical potential
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction
- K::Array{Float64,1}: Array of superlattice reciprocal wavevectors (noted K)

#Returns:

- Cumulant of the lattice Green's function (Complex{Float64}-valued matrix of size Nc x Nc)

"""
function periodize_cum_coex_suite(model::Model, kx::Float64, ky::Float64, K::Array{Float64,1})
    nambu_cum_inv = inv(periodize_AFM_orb_cum(model.cumulant_,kx,ky,K))
    k = [kx,ky]
    nambu_cum_inv[1,1]-=zetak(model,(k + K)...)+omegak(model,(k + K)...); nambu_cum_inv[2,2]-=zetak(model,(k + K)...)+omegak(model,(k + K)...); nambu_cum_inv[1,2]-=epsilonk(model,(k + K)...); nambu_cum_inv[2,1]-=epsilonk(model,(k + K)...)
    nambu_cum_inv[3,3]+=zetak(model,(k + K)...)+omegak(model,(k + K)...); nambu_cum_inv[4,4]+=zetak(model,(k + K)...)+omegak(model,(k + K)...); nambu_cum_inv[4,3]+=epsilonk(model,(k + K)...); nambu_cum_inv[3,4]+=epsilonk(model,(k + K)...)

    nambu_periodized_final = inv(nambu_cum_inv)
    return nambu_periodized_final
end

"""
Function associated to the function periodizing the cumulant in the case of coexistence between AF and SC. Last inversion step corresping to the
full Green's function.

#Arguments:

- model: model instance to have access to the attributes such as hopping elements, self-energy and chemical potential
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction

#Returns:

- nambu_periodized: Periodized lattice Green's function (Complex{Float64}-valued matrix Nc X Nc)

"""
function periodize_coex_cum_finale(model::Model, kx::Float64, ky::Float64)
    nambu_periodized = zeros(Complex{Float64},(4,4))
    Kx = [[0,0],[pi,0]]; Ky = [[0,0],[0,pi]]

    for K in Kx ## or Ky
        nambu_periodized += periodize_cum_coex_suite(model,kx,ky,K)
    end

    return 0.5*nambu_periodized
end

"""
Function to compute the superfluid stiffness in the case dSC and AF coexist (with cumulant periodization) perpendicular to the CuO2 cuprates plane
neglecting the current vertex corrections.

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction

#Returns:

- Superfluid stiffness: Value of the superfluid stiffness perpendicular to the plane for dSC+AF data with cumulant periodization

"""
function stiffness_cumulant_AFM_SC(model::Model, kx::Float64, ky::Float64)
    nambu_periodized_final = periodize_coex_cum_finale(model,kx,ky)
    return 1.0*real(tperp(kx,ky)*4.0*(nambu_periodized_final[2,4]*nambu_periodized_final[3,1]+nambu_periodized_final[1,4]*nambu_periodized_final[3,2]+
		         nambu_periodized_final[2,3]*nambu_periodized_final[4,1]+nambu_periodized_final[1,3]*nambu_periodized_final[4,2]))
end

"""
Function to compute the superfluid stiffness in the case there is only dSC (with cumulant periodization) perpendicular to the CuO2 cuprates plane
neglecting the current vertex corrections.

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction

#Returns:

- Superfluid stiffness: Value of the superfluid stiffness perpendicular to the plane for dSC-only data with cumulant periodization

"""
function stiffness_cum_SC(model::Model, kx::Float64, ky::Float64)
    nambu_periodized_final = periodize_nocoex_cum_finale(model,kx,ky)
    return 1.0*real(tperp(kx,ky)*4.0*(nambu_periodized_final[1,2]*nambu_periodized_final[2,1]))
end

"""
Function to compute the superfluid stiffness in the case there is only dSC (with Green periodization) perpendicular to the CuO2 cuprates plane and
neglecting the current vertex corrections

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction

#Returns:

- Superfluid stiffness: Value of the superfluid stiffness perpendicular to the ab plane for dSC-only data with Green periodization

"""
function stiffness_SC(model::Model, kx::Float64, ky::Float64)
    nambu_periodized = periodize_nocoex(kx, ky, build_gf_ktilde(model, kx, ky)) ###### NEW or not NEW <-------------
    return 1.0*real(4.0*tperp(kx,ky)*(nambu_periodized[1,2]*nambu_periodized[2,1])) ##Removed minus sign
end

"""
Function to compute the superfluid stiffness in the case there is only dSC (with Green periodization) in the CuO2 cuprates plane
neglecting the current vertex corrections

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction
- inplane_axis::String: String-valued argument to specify in-plane axis along which computing the superfluid stiffness (\'xx\' or \'yy\')
- AFM_SC_1::Int64: Integer value set to 1 to compute pure SC superfluid stiffness using the formula in coexisting AF+dSC regime. Otherwize, let it be 0

#Returns:

- Superfluid stiffness: Value of the superfluid stiffness in the ab plane for dSC-only data with Green periodization

"""
function stiffness_nocoex_per_ab(model::Model, kx::Float64, ky::Float64, inplane_axis::String, AFM_SC_1::Int64)
    if AFM_SC_1 == 1
        G = periodize_AFM_orb(build_gf_ktilde(model, kx, ky), kx, ky)
        if inplane_axis == "xx"
            kxky = [DxDxZetakbar,DxDxZetakbarepsilonk,DxDxepsilonk]
        elseif inplane_axis == "yy"
            kxky = [DyDyZetakbar,DyDyZetakbarepsilonk,DyDyepsilonk]
        else
            println("OUPS!") && throw("inplane_axis parameter takes only two possible values: \"xx\" or \"yy\".")
        end
        return 1.0*real(kxky[1](model,kx,ky)*4.0*(G[1,3]*G[3,1]+G[2,3]*G[3,2]+G[1,4]*G[4,1]+G[2,4]*G[4,2])+
                        kxky[2](model,kx,ky)*2.0*((G[1,3]+G[2,4])*(G[3,2]+G[4,1])+(G[1,4]+G[2,3])*(G[3,1]+G[4,2]))+
                        kxky[2](model,kx,ky)*2.0*((G[1,3]+G[2,4])*(G[3,2]+G[4,1])+(G[1,4]+G[2,3])*(G[3,1]+G[4,2]))+
                        kxky[3](model,kx,ky)*4.0*(G[2,4]*G[3,1]+G[1,4]*G[3,2]+G[2,3]*G[4,1]+G[1,3]*G[4,2]))
    elseif AFM_SC_1 == 0
        G = periodize_nocoex(kx, ky, build_gf_ktilde(model, kx, ky))
        if inplane_axis == "xx"
            kxky = DxDxEpsilonbark
        elseif inplane_axis == "yy"
            kxky = DyDyEpsilonbark
        else
            println("OUPS!") && throw("inplane_axis parameter takes only two possible values: \"xx\" or \"yy\".")
        end
        return 1.0*real(kxky(model,kx,ky)*4.0*(G[1,2]*G[2,1]))
    end
end

"""
Function to compute the superfluid stiffness in the case dSC+AF coexist (with Green periodization) perpendicular to the CuO2 cuprates plane and
neglecting the current vertex corrections

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction

#Returns:

- Superfluid stiffness: Value of the superfluid stiffness perpendicular to the ab plane for dSC+AF data with Green periodization

"""
function stiffness_test(model::Model, kx::Float64, ky::Float64)
    nambu_periodized = periodize_AFM_orb(build_gf_ktilde(model, kx, ky), kx, ky)
    return 1.0*real(4.0*tperp(kx,ky)*(nambu_periodized[2,4]*nambu_periodized[3,1]+nambu_periodized[1,4]*nambu_periodized[3,2]+
		           nambu_periodized[2,3]*nambu_periodized[4,1]+nambu_periodized[1,3]*nambu_periodized[4,2]))
end

"""
Function to compute the superfluid stiffness in the case dSC+AF coexist (with Green periodization) in the CuO2 cuprates plane
neglecting the current vertex corrections

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction
- inplane_axis::String: String-valued argument to specify in-plane axis along which computing the superfluid stiffness (\'xx\' or \'yy\')
- AFM_SC_1::Int64: Integer value set to 1 to compute pure SC superfluid stiffness using the formula in coexisting AF+dSC regime. Otherwise let it be 0

#Returns:

- Superfluid stiffness: Value of the superfluid stiffness in the ab plane for dSC+AF data with Green periodization

"""
function stiffness_coex_per_ab(model::Model, kx::Float64, ky::Float64, inplane_axis::String)
    G = periodize_AFM_orb(build_gf_ktilde(model, kx, ky), kx, ky)
    if inplane_axis == "xx"
        kxky = [DxDxEpsilonbark,DxDxZetakbar,DxDxZetakbarepsilonk,DxDxepsilonk]
    elseif inplane_axis == "yy"
        kxky = [DyDyEpsilonbark,DyDyZetakbar,DyDyZetakbarepsilonk,DyDyepsilonk]
    else
        println("OUPS!") && throw("inplane_axis parameter takes only two possible values: \"xx\" or \"yy\".")
    end
    return 1.0*real(kxky[1](model,kx,ky)*4.0*(G[1,3]*G[3,1]+G[2,3]*G[3,2]+G[1,4]*G[4,1]+G[2,4]*G[4,2])+
                    kxky[2](model,kx,ky)*2.0*((G[1,3]+G[2,4])*(G[3,2]+G[4,1])+(G[1,4]+G[2,3])*(G[3,1]+G[4,2]))+
                    kxky[2](model,kx,ky)*2.0*((G[1,3]+G[2,4])*(G[3,2]+G[4,1])+(G[1,4]+G[2,3])*(G[3,1]+G[4,2]))+
                    kxky[3](model,kx,ky)*4.0*(G[2,4]*G[3,1]+G[1,4]*G[3,2]+G[2,3]*G[4,1]+G[1,3]*G[4,2]))
end

"""
Function to compute the superfluid stiffness in the case there is dSC and AF (without periodization) perpendicular to the CuO2 cuprates plane
neglecting the current vertex corrections

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction

#Returns:

- Superfluid stiffness: Value of the superfluid stiffness perpendicular to the ab plane for dSC+AF data without periodization (tracing)

"""
function stiffness_cluster_G(model::Model, kx::Float64, ky::Float64)
    G = build_gf_ktilde(model,kx,ky)
    return 1.0*real(tperp(kx,ky)*4.0*(sum(G[1:4,5:end]))*(sum(G[5:end,1:4])))
end

"""
Function used to compute the superfluid stiffness in the case there is dSC and AF (without periodization) perpendicular to the CuO2 cuprates plane
neglecting the current vertex corrections

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction

#Returns:

- Superfluid stiffness: Value of the superfluid stiffness perpendicular to the ab plane for dSC+AF data without periodization

"""
function stiffness_cluster_G_first_neighbor(model::Model, kx::Float64, ky::Float64)
    G = build_gf_ktilde(model,kx,ky)
    return 1.0*real(tperp(kx,ky)*4.0*((G[3,7]+G[3,8]+G[4,7]+G[4,8])*(G[5,1]+G[5,2]+G[6,1]+G[6,2])+(G[1,7]+G[1,8]+
                    G[2,7]+G[2,8])*(G[5,3]+G[5,4]+G[6,3]+G[6,4])+(G[3,5]+G[3,6]+G[4,5]+G[4,6])*(G[7,1]+G[7,2]+G[8,1]+G[8,2])+(G[1,5]+G[1,6]+
                    G[2,5]+G[2,6])*(G[7,3]+G[7,4]+G[8,3]+G[8,4])))
end

"""
Function used to compute the superfluid stiffness in the case there is dSC and AF (without periodization) in the CuO2 cuprates plane
neglecting the current vertex corrections

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction
- inplane_axis::String: In-plane axis along which the superfluid stiffness is to be computed.

#Returns:

- Superfluid stiffness: Value of the superfluid stiffness in the ab plane for dSC+AF data without periodization

"""
function stiffness_cluster_G_ab(model::Model, kx::Float64, ky::Float64, inplane_axis::String)
    G = build_gf_ktilde(model,kx,ky)
    if inplane_axis == "xx"
        kxky = [DxDxomegak, DxDxzetakomegak, DxDxomegakepsilonk, DxDxzetakepsilonk, DxDxzetak, DxDxzetakepsilonk, DxDxomegakepsilonk, DxDxzetakepsilonk, DxDxepsilonk]
    elseif inplane_axis == "yy"
        kxky = [DyDyomegak, DyDyzetakomegak, DyDyomegakepsilonk, DyDyzetakepsilonk, DyDyzetak, DyDyzetakepsilonk, DyDyomegakepsilonk, DyDyzetakepsilonk, DyDyepsilonk]
    else
        println("OUPS!") && throw("inplane_axis parameter takes only two possible values: \"xx\" or \"yy\".")
    end
    return 1.0*real((kxky[1](model,kx,ky)*4.0*(G[1,5]*G[5,1]+G[2,5]*G[5,2]+G[3,5]*G[5,3]+G[4,5]*G[5,4]+G[1,6]*G[6,1]+G[2,6]*G[6,2]+G[3,6]*G[6,3]+G[4,6]*G[6,4]+G[1,7]*G[7,1]+G[2,7]*G[7,2]+G[3,7]*G[7,3]+G[4,7]*G[7,4]+G[1,8]*G[8,1]+G[2,8]*G[8,2]+G[3,8]*G[8,3]+G[4,8]*G[8,4])+
    kxky[2](model,kx,ky)*2.0*((G[1,5]+G[2,6])*(G[5,2]+G[6,1])+(G[1,6]+G[2,5])*(G[5,1]+G[6,2])+(G[3,5]+G[4,6])*(G[5,4]+G[6,3])+(G[3,6]+G[4,5])*(G[5,3]+G[6,4])+(G[1,7]+G[2,8])*(G[7,2]+G[8,1])+(G[1,8]+G[2,7])*(G[7,1]+G[8,2])+(G[3,7]+G[4,8])*(G[7,4]+G[8,3])+(G[3,8]+G[4,7])*(G[7,3]+G[8,4]))+
    kxky[3](model,kx,ky)*2.0*(G[3,5]*G[5,1]+G[4,5]*G[5,1]+G[3,5]*G[5,2]+G[4,5]*G[5,2]+G[3,7]*G[5,3]+G[3,8]*G[5,3]+(G[2,7]+G[2,8])*(G[5,2]+G[6,2])+(G[3,6]+G[4,6])*(G[6,1]+G[6,2])+G[3,7]*G[6,3]+G[3,8]*G[6,3]+(G[4,7]+G[4,8])*(G[5,4]+G[6,4])+G[2,6]*(G[6,3]+G[6,4])+G[3,7]*G[7,1]+G[4,7]*G[7,1]+G[3,7]*G[7,2]+G[4,7]*G[7,2]+G[3,5]*G[7,3]+G[3,6]*G[7,3]+G[4,5]*G[7,4]+G[4,6]*G[7,4]+G[2,7]*(G[7,3]+G[7,4])+G[1,7]*(G[5,1]+G[6,1]+G[7,3]+G[7,4])+G[3,8]*G[8,1]+G[4,8]*G[8,1]+G[1,5]*(G[5,3]+G[5,4]+G[7,1]+G[8,1])+G[1,6]*(G[6,3]+G[6,4]+G[7,1]+G[8,1])+(G[3,8]+G[4,8])*G[8,2]+G[2,6]*(G[7,2]+G[8,2])+G[2,5]*(G[5,3]+G[5,4]+G[7,2]+G[8,2])+G[3,5]*G[8,3]+G[3,6]*G[8,3]+(G[4,5]+G[4,6])*G[8,4]+G[2,8]*(G[8,3]+G[8,4])+G[1,8]*(G[5,1]+G[6,1]+G[8,3]+G[8,4]))+
    kxky[4](model,kx,ky)*2.0*((G[1,5]+G[2,6])*(G[5,2]+G[6,1])+(G[1,6]+G[2,5])*(G[5,1]+G[6,2])+(G[3,5]+G[4,6])*(G[5,4]+G[6,3])+(G[3,6]+G[4,5])*(G[5,3]+G[6,4])+(G[1,7]+G[2,8])*(G[7,2]+G[8,1])+(G[1,8]+G[2,7])*(G[7,1]+G[8,2])+(G[3,7]+G[4,8])*(G[7,4]+G[8,3])+(G[3,8]+G[4,7])*(G[7,3]+G[8,4]))+
    kxky[5](model,kx,ky)*4.0*(G[2,6]*G[5,1]+G[1,6]*G[5,2]+G[4,6]*G[5,3]+G[3,6]*G[5,4]+G[2,5]*G[6,1]+G[1,5]*G[6,2]+G[4,5]*G[6,3]+G[3,5]*G[6,4]+G[2,8]*G[7,1]+G[1,8]*G[7,2]+G[4,8]*G[7,3]+G[3,8]*G[7,4]+G[2,7]*G[8,1]+G[1,7]*G[8,2]+G[4,7]*G[8,3]+G[3,7]*G[8,4])+
    kxky[6](model,kx,ky)*2.0*(G[3,6]*G[5,1]+G[4,6]*G[5,1]+G[3,6]*G[5,2]+G[4,6]*G[5,2]+G[4,7]*G[5,3]+G[4,8]*G[5,3]+G[3,7]*G[5,4]+G[3,8]*G[5,4]+G[3,5]*G[6,1]+G[4,5]*G[6,1]+G[3,5]*G[6,2]+G[4,5]*G[6,2]+G[4,7]*G[6,3]+G[4,8]*G[6,3]+G[3,7]*G[6,4]+G[3,8]*G[6,4]+G[3,8]*G[7,1]+G[4,8]*G[7,1]+G[3,8]*G[7,2]+G[4,8]*G[7,2]+G[4,5]*G[7,3]+G[4,6]*G[7,3]+G[3,5]*G[7,4]+G[3,6]*G[7,4]+G[2,8]*(G[5,1]+G[6,1]+G[7,3]+G[7,4])+G[1,8]*(G[5,2]+G[6,2]+G[7,3]+G[7,4])+G[3,7]*G[8,1]+G[4,7]*G[8,1]+G[2,6]*(G[5,3]+G[5,4]+G[7,1]+G[8,1])+G[2,5]*(G[6,3]+G[6,4]+G[7,1]+G[8,1])+(G[3,7]+G[4,7])*G[8,2]+G[1,6]*(G[5,3]+G[5,4]+G[7,2]+G[8,2])+G[1,5]*(G[6,3]+G[6,4]+G[7,2]+G[8,2])+G[4,5]*G[8,3]+G[4,6]*G[8,3]+(G[3,5]+G[3,6])*G[8,4]+G[2,7]*(G[5,1]+G[6,1]+G[8,3]+G[8,4])+G[1,7]*(G[5,2]+G[6,2]+G[8,3]+G[8,4]))+
    kxky[7](model,kx,ky)*2.0*(G[3,5]*G[5,1]+G[4,5]*G[5,1]+G[3,5]*G[5,2]+G[4,5]*G[5,2]+G[3,7]*G[5,3]+G[3,8]*G[5,3]+(G[2,7]+G[2,8])*(G[5,2]+G[6,2])+(G[3,6]+G[4,6])*(G[6,1]+G[6,2])+G[3,7]*G[6,3]+G[3,8]*G[6,3]+(G[4,7]+G[4,8])*(G[5,4]+G[6,4])+G[2,6]*(G[6,3]+G[6,4])+G[3,7]*G[7,1]+G[4,7]*G[7,1]+G[3,7]*G[7,2]+G[4,7]*G[7,2]+G[3,5]*G[7,3]+G[3,6]*G[7,3]+G[4,5]*G[7,4]+G[4,6]*G[7,4]+G[2,7]*(G[7,3]+G[7,4])+G[1,7]*(G[5,1]+G[6,1]+G[7,3]+G[7,4])+G[3,8]*G[8,1]+G[4,8]*G[8,1]+G[1,5]*(G[5,3]+G[5,4]+G[7,1]+G[8,1])+G[1,6]*(G[6,3]+G[6,4]+G[7,1]+G[8,1])+(G[3,8]+G[4,8])*G[8,2]+G[2,6]*(G[7,2]+G[8,2])+G[2,5]*(G[5,3]+G[5,4]+G[7,2]+G[8,2])+G[3,5]*G[8,3]+G[3,6]*G[8,3]+(G[4,5]+G[4,6])*G[8,4]+G[2,8]*(G[8,3]+G[8,4])+G[1,8]*(G[5,1]+G[6,1]+G[8,3]+G[8,4]))+
    kxky[8](model,kx,ky)*2.0*(G[3,6]*G[5,1]+G[4,6]*G[5,1]+G[3,6]*G[5,2]+G[4,6]*G[5,2]+G[4,7]*G[5,3]+G[4,8]*G[5,3]+G[3,7]*G[5,4]+G[3,8]*G[5,4]+G[3,5]*G[6,1]+G[4,5]*G[6,1]+G[3,5]*G[6,2]+G[4,5]*G[6,2]+G[4,7]*G[6,3]+G[4,8]*G[6,3]+G[3,7]*G[6,4]+G[3,8]*G[6,4]+G[3,8]*G[7,1]+G[4,8]*G[7,1]+G[3,8]*G[7,2]+G[4,8]*G[7,2]+G[4,5]*G[7,3]+G[4,6]*G[7,3]+G[3,5]*G[7,4]+G[3,6]*G[7,4]+G[2,8]*(G[5,1]+G[6,1]+G[7,3]+G[7,4])+G[1,8]*(G[5,2]+G[6,2]+G[7,3]+G[7,4])+G[3,7]*G[8,1]+G[4,7]*G[8,1]+G[2,6]*(G[5,3]+G[5,4]+G[7,1]+G[8,1])+G[2,5]*(G[6,3]+G[6,4]+G[7,1]+G[8,1])+(G[3,7]+G[4,7])*G[8,2]+G[1,6]*(G[5,3]+G[5,4]+G[7,2]+G[8,2])+G[1,5]*(G[6,3]+G[6,4]+G[7,2]+G[8,2])+G[4,5]*G[8,3]+G[4,6]*G[8,3]+(G[3,5]+G[3,6])*G[8,4]+G[2,7]*(G[5,1]+G[6,1]+G[8,3]+G[8,4])+G[1,7]*(G[5,2]+G[6,2]+G[8,3]+G[8,4]))+
    kxky[9](model,kx,ky)*4.0*((G[3,7]+G[3,8]+G[4,7]+G[4,8])*(G[5,1]+G[5,2]+G[6,1]+G[6,2])+(G[1,7]+G[1,8]+G[2,7]+G[2,8])*(G[5,3]+G[5,4]+G[6,3]+G[6,4])+(G[3,5]+G[3,6]+G[4,5]+G[4,6])*(G[7,1]+G[7,2]+G[8,1]+G[8,2])+(G[1,5]+G[1,6]+G[2,5]+G[2,6])*(G[7,3]+G[7,4]+G[8,3]+G[8,4]))))
end

"""
Function to compute the superfluid stiffness in the case there is dSC only (without periodization) perpendicular to the CuO2 cuprates plane
neglecting the current vertex corrections

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction

#Returns:

- Superfluid stiffness: Value of the superfluid stiffness perpendicular to the ab plane for dSC-only data without periodization

"""
function stiffness_trace_G(model::Model, kx::Float64, ky::Float64)
    G = build_gf_ktilde(model,kx,ky)
    return 1.0*real(tperp(kx,ky)*4.0*trace(G[1:4,5:end]*G[5:end,1:4]))
end

"""
Function to compute the superfluid stiffness in the case there is dSC only (without periodization) in the CuO2 cuprates plane neglecting the current vertex corrections

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- kx::Float64: Wavevector in x direction
- ky::Float64: Wavevector in y direction

#Returns:

- Superfluid stiffness: Value of the superfluid stiffness in the ab plane for dSC-only data without periodization

"""
function stiffness_trace_G_ab(model::Model, kx::Float64, ky::Float64, inplane_axis::String)
    G = build_gf_ktilde(model,kx,ky)
    if inplane_axis == "xx"
        kxky = DxDxEpsilonbark
    elseif inplane_axis == "yy"
        kxky = DyDyEpsilonbark
    else
        println("OUPS!") && throw("inplane_axis parameter takes only two possible values: \"xx\" or \"yy\".")
    end
    return 1.0*real(kxky(model,kx,ky)*kxky(model,kx,ky)*4.0*trace(G[1:4,5:end]*G[5:end,1:4]))
end

"""
Template function useful to perform integration over the BZ when computing in-plane superfluid stiffness without periodization.
Called in module Stiffness.jl.

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- param::Tuple{String,String}: Tuple of two strings specifying the direction along which to compute the superfluid stiffness and the state (pure or mixed).

#Returns:

- Superfluid stiffness integrand. It is a function object used to compute integrals.

"""
function make_stiffness_cluster_G_kintegrand_ab(model::PeriodizeSC.Model, param::T) where T <: Tuple{AbstractString,AbstractString}
    function stiffness_kintegrand(kk::Array{Float64,1})
        if isa(param,Tuple{String,String})
            if param[1] == "xx" && param[2] == "COEX"
                return(PeriodizeSC.stiffness_cluster_G_ab(model,kk[1],kk[2],"xx"))
            elseif param[1] == "yy" && param[2] == "COEX"
                return(PeriodizeSC.stiffness_cluster_G_ab(model,kk[1],kk[2],"yy"))
            elseif param[1] == "xx" && param[2] == "NOCOEX"
                return(PeriodizeSC.stiffness_trace_G_ab(model,kk[1],kk[2],"xx"))
            elseif param[1] == "yy" && param[2] == "NOCOEX"
                return(PeriodizeSC.stiffness_trace_G_ab(model,kk[1],kk[2],"yy"))
            end
        else
            throw(ErrorException("Error occurred in function PeriodizeSC.make_stiffness_cluster_G_kintegrand_ab. Type not permitted in template."))
        end
    end
    return stiffness_kintegrand
end

"""
Template function useful to perform integration over the BZ when computing in-plane superfluid stiffness with periodization.
Called in module Stiffness.jl.

#Arguments:

- model: model instance to have access to the attributes such as hopping elements
- param::Tuple{String,String,Int64}: Tuple of three elements specifying the direction along which to compute the superfluid stiffness, the state (pure or mixed) and an internal debug parameter.

#Returns:

- Superfluid stiffness integrand. It is a function object used to compute integrals.

"""
function make_stiffness_kintegrand_per_ab(model::PeriodizeSC.Model, param::T) where T <: Tuple{AbstractString,AbstractString,Number}
    function stiffness_kintegrand(kk::Array{Float64,1})
        if isa(param,Tuple{String,String,Int64})
            if param[1] == "xx" && param[2] == "COEX"
                return(PeriodizeSC.stiffness_coex_per_ab(model,kk[1],kk[2],"xx"))
            elseif param[1] == "yy" && param[2] == "COEX"
                return(PeriodizeSC.stiffness_coex_per_ab(model,kk[1],kk[2],"yy"))
            elseif param[1] == "xx" && param[2] == "NOCOEX" && param[3] == 0
                return(PeriodizeSC.stiffness_nocoex_per_ab(model,kk[1],kk[2],"xx",0))
            elseif param[1] == "xx" && param[2] == "NOCOEX" && param[3] == 1
                return(PeriodizeSC.stiffness_nocoex_per_ab(model,kk[1],kk[2],"xx",1))
            elseif param[1] == "yy" && param[2] == "NOCOEX" && param[3] == 0
                return(PeriodizeSC.stiffness_nocoex_per_ab(model,kk[1],kk[2],"yy",0))
            elseif param[1] == "yy" && param[2] == "NOCOEX" && param[3] == 1
                return(PeriodizeSC.stiffness_nocoex_per_ab(model,kk[1],kk[2],"yy",1))
            end
        else
            throw(ErrorException("Error occurred in function PeriodizeSC.make_stiffness_cluster_G_kintegrand_ab. Type not permitted in template."))
        end
    end
    return stiffness_kintegrand
end

############################################################################## Important Note ###########################################################################
## The following functions act as decorators to the superfluid stiffness functions. These functions are necessary to compute the integrals over k-space. The situation ##
## in which each and every one of the functions is used is specified by its name.                                                                                ##
#########################################################################################################################################################################

function make_stiffness_trace_G_kintegrand(model::PeriodizeSC.Model)
    function stiffness_kintegrand(kk::Array{Float64,1})
        return(PeriodizeSC.stiffness_trace_G(model,kk[1],kk[2]))
    end
    return stiffness_kintegrand
end

function make_stiffness_cluster_G_kintegrand(model::PeriodizeSC.Model)
    function stiffness_kintegrand(kk::Array{Float64,1})
        return(PeriodizeSC.stiffness_cluster_G(model,kk[1],kk[2]))
    end
    return stiffness_kintegrand
end

function make_stiffness_cluster_G_kintegrand_fourth(model::PeriodizeSC.Model)
    function stiffness_kintegrand(kk::Array{Float64,1})
        return(PeriodizeSC.stiffness_cluster_G_first_neighbor(model,kk[1],kk[2]))
    end
    return stiffness_kintegrand
end

function make_stiffness_kintegrand_cum_AFM_SC(model::PeriodizeSC.Model)
    function stiffness_kintegrand(kk::Array{Float64, 1})
        return(PeriodizeSC.stiffness_cumulant_AFM_SC(model, kk[1], kk[2]))
    end
    return stiffness_kintegrand
end

function make_stiffness_kintegrand_cum_SC(model::PeriodizeSC.Model)
    function stiffness_kintegrand(kk::Array{Float64, 1})
        return(PeriodizeSC.stiffness_cum_SC(model, kk[1], kk[2]))
    end
    return stiffness_kintegrand
end

function make_stiffness_kintegrand_SC(model::PeriodizeSC.Model)
    function stiffness_kintegrand(kk::Array{Float64,1})
        return(PeriodizeSC.stiffness_SC(model, kk[1], kk[2]))
    end
    return stiffness_kintegrand
end

function make_stiffness_kintegrand_test(model::PeriodizeSC.Model)
    function stiffness_kintegrand(kk::Array{Float64,1})
        return(PeriodizeSC.stiffness_test(model,kk[1],kk[2]))
    end
    return stiffness_kintegrand
end

###################################################################### End of Important Note ################################################################################
#############################################################################################################################################################################

"""
Function used to sum over a grid of the reduced BZ.

#Arguments:

- modelvector: modelvector instance
- fct: Function object to be used to sum over the reduced BZ.
- gridK::Int64: k-space grid NxN, where N is set to 100 by default.
- len_sEvec_c::Int64: Length of the Matsubara frequency grid. Set by default to 500.

#Returns:

- Returns an array of two columns contaning the k-integrated Green's function and the associated Matsubara frequency.

"""
function sum_RBZ(modelvector::ModelVector, fct; gridK::Int64 = 100, len_sEvec_c::Int64=500)
    len_sEvec_c > size(modelvector.sEvec_c_)[1] && throw(ErrorException("You have just exceeded the number of Matsubara frequencies available for calculations!"))
    kx = linspace(-pi,pi,gridK)
    ky = linspace(-pi,pi,gridK)
    println("Length of len_sEvec_c and length of modelvector.sEvec_c_, respectively :", len_sEvec_c, size(modelvector.sEvec_c_)[1])
    result_k_ind = Array{Float64}(len_sEvec_c)
    for n in 1:len_sEvec_c
        stiffness = 0.0
        model = Model(modelvector,n)
        for i in 1:size(ky)[1]
            for j in 1:size(kx)[1]
                stiffness += fct(model,kx[j],ky[i])
            end
        end
        stiffness = 1./(2.0*gridK^2)*stiffness
        result_k_ind[n] = stiffness
        println("stiffness = ", stiffness)
    end
    result_out = hcat(modelvector.wvec_[1:len_sEvec_c], result_k_ind)
    return(result_out)
end

function sum_BZ(modelvector::ModelVector, fct; gridK::Int64 = 100, len_sEvec_c::Int64=500)
    println("in sum_BZ: 2.0*sum_RBZ")
    result_out = 2.0*sum_RBZ(modelvector, fct, gridK=gridK, len_sEvec_c=len_sEvec_c)
    return result_out
end

"""
Function used to integrate over the BZ when computing out-of-plane superfluid stiffness.

#Arguments:Integer value set to 1 to compute pure SC superfluid stiffness using the formula in coexisting AF+dSC regime. Otherwize, let it be 0

- modelvectInteger value set to 1 to compute pure SC superfluid stiffness using the formula in coexisting AF+dSC regime. Otherwize, let it be 0
- fct: FuncInteger value set to 1 to compute pure SC superfluid stiffness using the formula in coexisting AF+dSC regime. Otherwize, let it be 0
- gridK::InInteger value set to 1 to compute pure SC superfluid stiffness using the formula in coexisting AF+dSC regime. Otherwize, let it be 0
- len_sEvecInteger value set to 1 to compute pure SC superfluid stiffness using the formula in coexisting AF+dSC regime. Otherwize, let it be 0

#Returns:

- Returns an array of two columns contaning the k-integrated Green's function and the associated Matsubara frequency.

"""
function calcintegral_RBZ(modelvector::ModelVector, fct; maxevals::Int64=100_000, len_sEvec_c::Int64=500)
    len_sEvec_c > size(modelvector.sEvec_c_)[1] && throw(ErrorException("You have just exceeded the number of Matsubara frequencies available for calculations!"))
    println("Length of len_sEvec_c and length of modelvector.sEvec_c_, respectively : ", len_sEvec_c, "  ", size(modelvector.sEvec_c_)[1])
    result = Array{Float64,1}(len_sEvec_c)
    println("in calcintegral_RBZ, kwargs")
    for n in 1:len_sEvec_c
        model = Model(modelvector, n)
        #println("w_: ", model.w_)
        result[n] = (1.0/2.0)*(2.0*pi)^(-2.0)*hcubature(fct(model), (-pi,-pi), (pi,pi), reltol=1.49e-8, abstol=1.49e-8, maxevals=maxevals)[1]
        #println("stiff: ", result[n])
    end
    result_out = hcat(modelvector.wvec_[1:len_sEvec_c], result)
    return result_out
end

function calcintegral_BZ(modelvector::ModelVector, fct; maxevals::Int64=100_000, len_sEvec_c::Int64=500)
    println("in calcintegral_BZ: 2.0*calcintegral_RBZ")
    result_out = 2.0*calcintegral_RBZ(modelvector,fct,len_sEvec_c=len_sEvec_c)
    return result_out
end

"""
Function used to integrate over the BZ when computing in-plane superfluid stiffness.

#Arguments:

- modelvector: modelvector instance
- fct: Function object to be used to integrate over the BZ.
- param::Tuple{String,String,Int64}: Tuple of three elements necessary to provide because of the form of the parameter fct used.
- gridK::Int64: k-space grid NxN, where N is set to 100 by default.
- len_sEvec_c::Int64: Length of the Matsubara frequency grid. Set by default to 500.

#Returns:

- Returns an array of two columns contaning the k-integrated Green's function and the associated Matsubara frequency.

"""
function calcintegral_RBZ_ab(modelvector::ModelVector, fct, param::T; maxevals::Int64=100_000, len_sEvec_c::Int64=500) where T <: Tuple
    len_sEvec_c > size(modelvector.sEvec_c_)[1] && throw(ErrorException("You have just exceeded the number of Matsubara frequencies available for calculations!"))
    println("Length of len_sEvec_c and length of modelvector.sEvec_c_, respectively : ", len_sEvec_c, "  ",size(modelvector.sEvec_c_)[1])
    result = Array{Float64,1}(len_sEvec_c)
    println("in calcintegral_RBZ_ab, kwargs")
    for n in 1:len_sEvec_c
        model = Model(modelvector, n)
        result[n] = (1.0/2.0)*(2.0*pi)^(-2.0)*hcubature(fct(model,param), (-pi,-pi), (pi,pi), reltol=1.49e-8, abstol=1.49e-8, maxevals=maxevals)[1]
        println("stiff: ", result[n])
    end
    result_out = hcat(modelvector.wvec_[1:len_sEvec_c], result)
    return result_out
end

function calcintegral_BZ_ab(modelvector::ModelVector, fct, param::T; maxevals::Int64=100_000, len_sEvec_c::Int64=500) where T <: Tuple
    println("in calcintegral_BZ_ab: 2.0*calcintegral_RBZ_ab")
    result_out = 2.0*calcintegral_RBZ_ab(modelvector,fct,param,maxevals=maxevals,len_sEvec_c=len_sEvec_c)
    return result_out
end

"""
Function used to map any functions on desired k-space grid. Used when pre-computing current vertices.

#Arguments:

- model: model instance
- grid::Int64: k-space grid dimension.
- fct: Function object to be mapped onto the BZ grid.

#Returns:

- Returns an array object of type Array{Float64}(Grid,Grid). This object contains the fct values on the BZ grid.

"""
function k_grid(model::Model,Grid::Int64,fct) #This function can take tktilde or dispersion relation
    try
        BZ_grid = Array{Float64}(Grid,Grid)
        for (i,ky) in enumerate(linspace(-pi,pi,Grid))
            for (j,kx) in enumerate(linspace(-pi,pi,Grid))
                BZ_grid[i,j] = fct(model,kx,ky)
            end
        end
        return BZ_grid
    catch ex
        if isa(ex,MethodError)
            BZ_grid = Array{Array{Complex{Float64},2}}(Grid,Grid)
            for (i,ky) in enumerate(linspace(-pi,pi,Grid))
                for (j,kx) in enumerate(linspace(-pi,pi,Grid))
                    BZ_grid[i,j] = fct(model,kx,ky)
                end
            end
            return BZ_grid
        else
            throw(ErrorException("Problem in definition of the type of BZ_grid. See function PeriodizeSC.k_grid."))
        end
    end
end

"""
Function used to in-plane AF+dSC superfluid stiffness without periodization.

#Arguments:

- modelvector: modelvector instance
- kxky::Array{Array{Float64,2},1}: Pre-computed current vertices.
- tk::Array{Array{Complex{Float64},2},2}: Pre-computed hopping matrix represented in the reduced supercluster BZ.
- Grid::Int64: Value of the k-space grid dimension.
- len_sEvec_c::Int64: Length of the Matsubara frequency grid. Set by default to 500.

#Returns:

- Returns an array of two columns contaning the k-integrated Green's function and the associated Matsubara frequency.

"""
function stiffness_cluster_G_ab_k_grid(modelvector::ModelVector, kxky::Array{Array{Float64,2},1}, tk::Array{Array{Complex{Float64},2},2}, Grid::Int64; len_sEvec_c::Int64=500)
    len_sEvec_c > size(modelvector.sEvec_c_)[1] && throw(ErrorException("You have exceeded the number of Matsubara frequencies available for calculations."))
    result_n = Array{Float64,1}(len_sEvec_c)
    for n in 1:len_sEvec_c
        model = Model(modelvector,n) ########################### Should devide by two when tracing ???
        Sum = 0.0
        for (i,ky) in enumerate(linspace(-pi,pi,Grid))
            for (j,kx) in enumerate(linspace(-pi,pi,Grid))
                G = build_gf_ktilde_prebuild(model,tk[i,j])
                Sum += 1.0*(kxky[1][i,j]*4.0*(G[1,5]*G[5,1]+G[2,5]*G[5,2]+G[3,5]*G[5,3]+G[4,5]*G[5,4]+G[1,6]*G[6,1]+G[2,6]*G[6,2]+G[3,6]*G[6,3]+G[4,6]*G[6,4]+G[1,7]*G[7,1]+G[2,7]*G[7,2]+G[3,7]*G[7,3]+G[4,7]*G[7,4]+G[1,8]*G[8,1]+G[2,8]*G[8,2]+G[3,8]*G[8,3]+G[4,8]*G[8,4])+
                kxky[2][i,j]*2.0*((G[1,5]+G[2,6])*(G[5,2]+G[6,1])+(G[1,6]+G[2,5])*(G[5,1]+G[6,2])+(G[3,5]+G[4,6])*(G[5,4]+G[6,3])+(G[3,6]+G[4,5])*(G[5,3]+G[6,4])+(G[1,7]+G[2,8])*(G[7,2]+G[8,1])+(G[1,8]+G[2,7])*(G[7,1]+G[8,2])+(G[3,7]+G[4,8])*(G[7,4]+G[8,3])+(G[3,8]+G[4,7])*(G[7,3]+G[8,4]))+
                kxky[3][i,j]*2.0*(G[3,5]*G[5,1]+G[4,5]*G[5,1]+G[3,5]*G[5,2]+G[4,5]*G[5,2]+G[3,7]*G[5,3]+G[3,8]*G[5,3]+(G[2,7]+G[2,8])*(G[5,2]+G[6,2])+(G[3,6]+G[4,6])*(G[6,1]+G[6,2])+G[3,7]*G[6,3]+G[3,8]*G[6,3]+(G[4,7]+G[4,8])*(G[5,4]+G[6,4])+G[2,6]*(G[6,3]+G[6,4])+G[3,7]*G[7,1]+G[4,7]*G[7,1]+G[3,7]*G[7,2]+G[4,7]*G[7,2]+G[3,5]*G[7,3]+G[3,6]*G[7,3]+G[4,5]*G[7,4]+G[4,6]*G[7,4]+G[2,7]*(G[7,3]+G[7,4])+G[1,7]*(G[5,1]+G[6,1]+G[7,3]+G[7,4])+G[3,8]*G[8,1]+G[4,8]*G[8,1]+G[1,5]*(G[5,3]+G[5,4]+G[7,1]+G[8,1])+G[1,6]*(G[6,3]+G[6,4]+G[7,1]+G[8,1])+(G[3,8]+G[4,8])*G[8,2]+G[2,6]*(G[7,2]+G[8,2])+G[2,5]*(G[5,3]+G[5,4]+G[7,2]+G[8,2])+G[3,5]*G[8,3]+G[3,6]*G[8,3]+(G[4,5]+G[4,6])*G[8,4]+G[2,8]*(G[8,3]+G[8,4])+G[1,8]*(G[5,1]+G[6,1]+G[8,3]+G[8,4]))+
                kxky[4][i,j]*2.0*((G[1,5]+G[2,6])*(G[5,2]+G[6,1])+(G[1,6]+G[2,5])*(G[5,1]+G[6,2])+(G[3,5]+G[4,6])*(G[5,4]+G[6,3])+(G[3,6]+G[4,5])*(G[5,3]+G[6,4])+(G[1,7]+G[2,8])*(G[7,2]+G[8,1])+(G[1,8]+G[2,7])*(G[7,1]+G[8,2])+(G[3,7]+G[4,8])*(G[7,4]+G[8,3])+(G[3,8]+G[4,7])*(G[7,3]+G[8,4]))+
                kxky[5][i,j]*4.0*(G[2,6]*G[5,1]+G[1,6]*G[5,2]+G[4,6]*G[5,3]+G[3,6]*G[5,4]+G[2,5]*G[6,1]+G[1,5]*G[6,2]+G[4,5]*G[6,3]+G[3,5]*G[6,4]+G[2,8]*G[7,1]+G[1,8]*G[7,2]+G[4,8]*G[7,3]+G[3,8]*G[7,4]+G[2,7]*G[8,1]+G[1,7]*G[8,2]+G[4,7]*G[8,3]+G[3,7]*G[8,4])+
                kxky[6][i,j]*2.0*(G[3,6]*G[5,1]+G[4,6]*G[5,1]+G[3,6]*G[5,2]+G[4,6]*G[5,2]+G[4,7]*G[5,3]+G[4,8]*G[5,3]+G[3,7]*G[5,4]+G[3,8]*G[5,4]+G[3,5]*G[6,1]+G[4,5]*G[6,1]+G[3,5]*G[6,2]+G[4,5]*G[6,2]+G[4,7]*G[6,3]+G[4,8]*G[6,3]+G[3,7]*G[6,4]+G[3,8]*G[6,4]+G[3,8]*G[7,1]+G[4,8]*G[7,1]+G[3,8]*G[7,2]+G[4,8]*G[7,2]+G[4,5]*G[7,3]+G[4,6]*G[7,3]+G[3,5]*G[7,4]+G[3,6]*G[7,4]+G[2,8]*(G[5,1]+G[6,1]+G[7,3]+G[7,4])+G[1,8]*(G[5,2]+G[6,2]+G[7,3]+G[7,4])+G[3,7]*G[8,1]+G[4,7]*G[8,1]+G[2,6]*(G[5,3]+G[5,4]+G[7,1]+G[8,1])+G[2,5]*(G[6,3]+G[6,4]+G[7,1]+G[8,1])+(G[3,7]+G[4,7])*G[8,2]+G[1,6]*(G[5,3]+G[5,4]+G[7,2]+G[8,2])+G[1,5]*(G[6,3]+G[6,4]+G[7,2]+G[8,2])+G[4,5]*G[8,3]+G[4,6]*G[8,3]+(G[3,5]+G[3,6])*G[8,4]+G[2,7]*(G[5,1]+G[6,1]+G[8,3]+G[8,4])+G[1,7]*(G[5,2]+G[6,2]+G[8,3]+G[8,4]))+
                kxky[7][i,j]*2.0*(G[3,5]*G[5,1]+G[4,5]*G[5,1]+G[3,5]*G[5,2]+G[4,5]*G[5,2]+G[3,7]*G[5,3]+G[3,8]*G[5,3]+(G[2,7]+G[2,8])*(G[5,2]+G[6,2])+(G[3,6]+G[4,6])*(G[6,1]+G[6,2])+G[3,7]*G[6,3]+G[3,8]*G[6,3]+(G[4,7]+G[4,8])*(G[5,4]+G[6,4])+G[2,6]*(G[6,3]+G[6,4])+G[3,7]*G[7,1]+G[4,7]*G[7,1]+G[3,7]*G[7,2]+G[4,7]*G[7,2]+G[3,5]*G[7,3]+G[3,6]*G[7,3]+G[4,5]*G[7,4]+G[4,6]*G[7,4]+G[2,7]*(G[7,3]+G[7,4])+G[1,7]*(G[5,1]+G[6,1]+G[7,3]+G[7,4])+G[3,8]*G[8,1]+G[4,8]*G[8,1]+G[1,5]*(G[5,3]+G[5,4]+G[7,1]+G[8,1])+G[1,6]*(G[6,3]+G[6,4]+G[7,1]+G[8,1])+(G[3,8]+G[4,8])*G[8,2]+G[2,6]*(G[7,2]+G[8,2])+G[2,5]*(G[5,3]+G[5,4]+G[7,2]+G[8,2])+G[3,5]*G[8,3]+G[3,6]*G[8,3]+(G[4,5]+G[4,6])*G[8,4]+G[2,8]*(G[8,3]+G[8,4])+G[1,8]*(G[5,1]+G[6,1]+G[8,3]+G[8,4]))+
                kxky[8][i,j]*2.0*(G[3,6]*G[5,1]+G[4,6]*G[5,1]+G[3,6]*G[5,2]+G[4,6]*G[5,2]+G[4,7]*G[5,3]+G[4,8]*G[5,3]+G[3,7]*G[5,4]+G[3,8]*G[5,4]+G[3,5]*G[6,1]+G[4,5]*G[6,1]+G[3,5]*G[6,2]+G[4,5]*G[6,2]+G[4,7]*G[6,3]+G[4,8]*G[6,3]+G[3,7]*G[6,4]+G[3,8]*G[6,4]+G[3,8]*G[7,1]+G[4,8]*G[7,1]+G[3,8]*G[7,2]+G[4,8]*G[7,2]+G[4,5]*G[7,3]+G[4,6]*G[7,3]+G[3,5]*G[7,4]+G[3,6]*G[7,4]+G[2,8]*(G[5,1]+G[6,1]+G[7,3]+G[7,4])+G[1,8]*(G[5,2]+G[6,2]+G[7,3]+G[7,4])+G[3,7]*G[8,1]+G[4,7]*G[8,1]+G[2,6]*(G[5,3]+G[5,4]+G[7,1]+G[8,1])+G[2,5]*(G[6,3]+G[6,4]+G[7,1]+G[8,1])+(G[3,7]+G[4,7])*G[8,2]+G[1,6]*(G[5,3]+G[5,4]+G[7,2]+G[8,2])+G[1,5]*(G[6,3]+G[6,4]+G[7,2]+G[8,2])+G[4,5]*G[8,3]+G[4,6]*G[8,3]+(G[3,5]+G[3,6])*G[8,4]+G[2,7]*(G[5,1]+G[6,1]+G[8,3]+G[8,4])+G[1,7]*(G[5,2]+G[6,2]+G[8,3]+G[8,4]))+
                kxky[9][i,j]*4.0*((G[3,7]+G[3,8]+G[4,7]+G[4,8])*(G[5,1]+G[5,2]+G[6,1]+G[6,2])+(G[1,7]+G[1,8]+G[2,7]+G[2,8])*(G[5,3]+G[5,4]+G[6,3]+G[6,4])+(G[3,5]+G[3,6]+G[4,5]+G[4,6])*(G[7,1]+G[7,2]+G[8,1]+G[8,2])+(G[1,5]+G[1,6]+G[2,5]+G[2,6])*(G[7,3]+G[7,4]+G[8,3]+G[8,4])))
            end
        end
        result_n[n] = 2./(Grid^2)*real(Sum) # Factor 2 accounts for the imaginary axis
    end
    result_out = hcat(modelvector.wvec_[1:len_sEvec_c],result_n)
    println(result_out)
    return result_out
end

"""
Function used to in-plane AF+dSC superfluid stiffness with periodization.

#Arguments:

- modelvector: modelvector instance
- kxky::Array{Array{Float64,2},1}: Pre-computed current vertices.
- tk::Array{Array{Complex{Float64},2},2}: Pre-computed hopping matrix represented in the reduced supercluster BZ.
- Grid::Int64: Value of the k-space grid dimension.
- Cum::Int64: 1 if the periodization is to be done on the cumulant. 0 if the periodization is to be done on Green's function.
- super_data_M_el::Float64: Value of the AF order parameter amplitude.
- M_tol::Float64: Relevant for debugging purposes only.
- len_sEvec_c::Int64: Length of the Matsubara frequency grid. Set by default to 500.

#Returns:

- Returns an array of two columns contaning the k-integrated Green's function and the associated Matsubara frequency.

"""
function stiffness_COEX_Per_Cum_ab_k_grid(modelvector::ModelVector, kxky::Array{Array{Float64,2},1}, tk::Array{Array{Complex{Float64},2},2}, Grid::Int64, Cum::Int64, super_data_M_el::Float64, M_tol::Float64; len_sEvec_c::Int64=500)
    len_sEvec_c > size(modelvector.sEvec_c_)[1] && throw(ErrorException("You have exceeded the number of Matsubara frequencies available for calculations."))
    result_n = Array{Float64,1}(len_sEvec_c)
    for n in 1:len_sEvec_c
        model = Model(modelvector,n)
        Sum = 0.0
        cond = abs(super_data_M_el) > M_tol
        if cond
            for (i,ky) in enumerate(linspace(-pi,pi,Grid))
                for (j,kx) in enumerate(linspace(-pi,pi,Grid))
                    if Cum == 0
                        gf_ktilde = build_gf_ktilde_prebuild(model,tk[i,j])
                        G = periodize_AFM_orb(gf_ktilde, kx, ky)
                    elseif Cum == 1
                        G = periodize_coex_cum_finale(model, kx, ky)
                    end
                    Sum += 1.0*(kxky[2][i,j]*4.0*(G[1,3]*G[3,1]+G[2,3]*G[3,2]+G[1,4]*G[4,1]+G[2,4]*G[4,2])+
                    kxky[3][i,j]*2.0*((G[1,3]+G[2,4])*(G[3,2]+G[4,1])+(G[1,4]+G[2,3])*(G[3,1]+G[4,2]))+
                    kxky[3][i,j]*2.0*((G[1,3]+G[2,4])*(G[3,2]+G[4,1])+(G[1,4]+G[2,3])*(G[3,1]+G[4,2]))+
                    kxky[4][i,j]*4.0*(G[2,4]*G[3,1]+G[1,4]*G[3,2]+G[2,3]*G[4,1]+G[1,3]*G[4,2]))
                end
            end
        else
            for (i,ky) in enumerate(linspace(-pi,pi,Grid))
                for (j,kx) in enumerate(linspace(-pi,pi,Grid))
                    if Cum == 0
                        gf_ktilde = build_gf_ktilde_prebuild(model,tk[i,j])
                        G = periodize_nocoex(kx, ky, gf_ktilde)
                    elseif Cum == 1
                        G = stiffness_cum_SC(model, kx, ky)
                    end
                    Sum += 2.0*(kxky[1][i,j]*4.0*(G[1,2]*G[2,1])) ## Factor 2.0 for the spin degrees of freedom
                end
            end
        end
        result_n[n] = (1./(2*Grid^2))*real(Sum) # Factor 1/2 to avoid counting the spin degree of freedom twice
    end
    result_out = hcat(modelvector.wvec_[1:len_sEvec_c],result_n)
    println(result_out)
    return result_out
end

"""
Function used to in-plane dSC superfluid stiffness with periodization.

#Arguments:

- modelvector: modelvector instance
- kxky::Array{Array{Float64,2},1}: Pre-computed current vertices.
- tk::Array{Array{Complex{Float64},2},2}: Pre-computed hopping matrix represented in the reduced supercluster BZ.
- Grid::Int64: Value of the k-space grid dimension.
- Cum::Int64: 1 if the periodization is to be done on the cumulant. 0 if the periodization is to be done on Green's function.
- AFMSC::Int64: Integer value set to 1 to compute pure SC superfluid stiffness using the formula in coexisting AF+dSC regime. Otherwize, let it be 0.
- len_sEvec_c::Int64: Length of the Matsubara frequency grid. Set by default to 500.

#Returns:

- Returns an array of two columns contaning the k-integrated Green's function and the associated Matsubara frequency.

"""
function stiffness_NOCOEX_Per_Cum_ab_k_grid(modelvector::ModelVector, kxky::Array{Array{Float64,2},1}, tk::Array{Array{Complex{Float64},2},2}, Grid::Int64, Cum::Int64, AFMSC::Int64; len_sEvec_c::Int64=500)
    len_sEvec_c > size(modelvector.sEvec_c_)[1] && throw(ErrorException("You have exceeded the number of Matsubara frequencies available for calculations."))
    result_n = Array{Float64,1}(len_sEvec_c)
    for n in 1:len_sEvec_c
        model = Model(modelvector,n)
        Sum = 0.0
        for (i,ky) in enumerate(linspace(-pi,pi,Grid))
            for (j,kx) in enumerate(linspace(-pi,pi,Grid))
                if Cum == 0 && AFMSC == 0
                    gf_ktilde = build_gf_ktilde_prebuild(model,tk[i,j])
                    G = periodize_nocoex(kx, ky, gf_ktilde)
                    Sum += 1.0*(kxky[1][i,j]*4.0*(G[1,2]*G[2,1]))
                elseif Cum == 1 && AFMSC == 0
                    G = periodize_cumulant_SC(model, kx, ky)
                    Sum += 1.0*(kxky[1][i,j]*4.0*(G[1,2]*G[2,1]))
                elseif Cum == 0 && AFMSC == 1
                    gf_ktilde = build_gf_ktilde_prebuild(model,tk[i,j])
                    G = periodize_AFM_orb(gf_ktilde, kx, ky)
                    Sum += 1.0*(kxky[1][i,j]*4.0*(G[1,3]*G[3,1]+G[2,3]*G[3,2]+G[1,4]*G[4,1]+G[2,4]*G[4,2])+
                    kxky[2][i,j]*2.0*((G[1,3]+G[2,4])*(G[3,2]+G[4,1])+(G[1,4]+G[2,3])*(G[3,1]+G[4,2]))+
                    kxky[2][i,j]*2.0*((G[1,3]+G[2,4])*(G[3,2]+G[4,1])+(G[1,4]+G[2,3])*(G[3,1]+G[4,2]))+
                    kxky[3][i,j]*4.0*(G[2,4]*G[3,1]+G[1,4]*G[3,2]+G[2,3]*G[4,1]+G[1,3]*G[4,2]))
                elseif Cum == 1 && AFMSC == 1
                    G = periodize_coex_cum_finale(model, kx, ky)
                    Sum += 1.0*(kxky[1][i,j]*4.0*(G[1,3]*G[3,1]+G[2,3]*G[3,2]+G[1,4]*G[4,1]+G[2,4]*G[4,2])+
                    kxky[2][i,j]*2.0*((G[1,3]+G[2,4])*(G[3,2]+G[4,1])+(G[1,4]+G[2,3])*(G[3,1]+G[4,2]))+
                    kxky[2][i,j]*2.0*((G[1,3]+G[2,4])*(G[3,2]+G[4,1])+(G[1,4]+G[2,3])*(G[3,1]+G[4,2]))+
                    kxky[3][i,j]*4.0*(G[2,4]*G[3,1]+G[1,4]*G[3,2]+G[2,3]*G[4,1]+G[1,3]*G[4,2]))
                end
            end
        end
        result_n[n] = 1./(Grid^2)*real(Sum) # No factor 2 to count the spin degree of freedom
    end
    result_out = hcat(modelvector.wvec_[1:len_sEvec_c],result_n)
    println(result_out)
    return result_out
end

end ## End of module

max_stars_repo_name: oliviersimard/SuperStiff

max_issues_repo_name: oliviersimard/SuperStiff

max_forks_repo_name: oliviersimard/SuperStiff


text:
@inline score( one::A, two::A  ) where {A <: UngappedAlignment} = @fastmath score( one ) + score( two )
@inline identity( one::A, two::A, onelen::Int, twolen::Int ) where {A <: UngappedAlignment} = @fastmath score( one, two ) / ( onelen + twolen )

@inline function sorted_offsets( sa::UnitRange, index::FMIndexes.FMIndex )
   res = Vector{Int}(length(sa))
   ridx = 1
   for sidx in sa
      res[ridx] = FMIndexes.sa_value( sidx, index ) + 1
      ridx += 1
   end
   sort!(res)
   res
end

# Paired version.
function splice_by_score!( one::Vector{A}, two::Vector{A}, threshold, buffer ) where A <: UngappedAlignment
   i = 1
   while i <= length( one )
      if threshold - score( one[i], two[i] ) > buffer
         splice!( one, i )
         splice!( two, i )
         i -= 1
      end
      i += 1
   end
end

function splice_by_identity!( one::Vector{A}, two::Vector{A}, threshold, buffer, onelen, twolen ) where A <: UngappedAlignment
   i = 1
   while i <= length( one )
      if threshold - identity( one[i], two[i], onelen, twolen ) > buffer
         splice!( one, i )
         splice!( two, i )
         i -= 1
      end
      i += 1
   end
end

# Paired ungapped_align
function ungapped_align( p::AlignParam, lib::GraphLib, fwd::FASTQRecord, rev::FASTQRecord; ispos=true, anchor_left=true )

   fwd_seed,fwd_readloc,strand = seed_locate( p, lib.index, fwd, offset_left=anchor_left, both_strands=!p.is_stranded )

   if (p.is_pair_rc && strand) || (!p.is_pair_rc && !strand)
      reverse_complement!(rev)
      rev_seed,rev_readloc = seed_locate( p, lib.index, rev, offset_left=!anchor_left, both_strands=false)
   else
      rev_seed,rev_readloc = seed_locate( p, lib.index, rev, offset_left=anchor_left, both_strands=false)
   end
   fwd_len = length(fwd.sequence)
   rev_len = length(rev.sequence)

   if !strand
      reverse_complement!(fwd)
      ispos = false
   end

   fwd_res      = Nullable{Vector{SGAlignment}}()
   rev_res      = Nullable{Vector{SGAlignment}}()
   maxscore = 0.0

   fwd_sorted = sorted_offsets( fwd_seed, lib.index )
   rev_sorted = sorted_offsets( rev_seed, lib.index )

   fidx,ridx = 1,1

   while fidx <= length(fwd_sorted) && ridx <= length(rev_sorted)

      dist = abs( fwd_sorted[fidx] - rev_sorted[ridx] )
      if dist < p.seed_pair_range # inside allowable pair distance

         geneind = convert( CoordInt, searchsortedlast( lib.offset, fwd_sorted[fidx] ) )
         next_offset= geneind < length(lib.offset) ? lib.offset[geneind+1] :
                                                           lib.offset[geneind] + length(lib.graphs[geneind].seq)
         if lib.offset[geneind] <= rev_sorted[ridx] < next_offset

            fwd_aln = _ungapped_align( p, lib, fwd, fwd_sorted[fidx], fwd_readloc; ispos=ispos, geneind=geneind )
            rev_aln = _ungapped_align( p, lib, rev, rev_sorted[ridx], rev_readloc; ispos=!ispos, geneind=geneind )

            @fastmath const scvar = identity( fwd_aln, rev_aln, fwd_len, rev_len )

            if isvalid(fwd_aln) && isvalid(rev_aln)
               if isnull( fwd_res ) || isnull( rev_res )
                  fwd_res = Nullable(SGAlignment[ fwd_aln ])
                  rev_res = Nullable(SGAlignment[ rev_aln ])
                  maxscore = scvar
               else
               # new best score
                  if scvar > maxscore
                     # better than threshold for all of the previous scores
                     if scvar - maxscore > p.score_range
                        if length( fwd_res.value ) >= 1
                           empty!( fwd_res.value )
                           empty!( rev_res.value )
                        end
                     else # keep at least one previous score
                        splice_by_identity!( fwd_res.value, rev_res.value, scvar, p.score_range, fwd_len, rev_len )
                     end
                     maxscore = scvar
                     push!( fwd_res.value, fwd_aln )
                     push!( rev_res.value, rev_aln )
                  else # new score is lower than previously seen
                     if maxscore - scvar <= p.score_range # but tolerable
                        push!( fwd_res.value, fwd_aln )
                        push!( rev_res.value, rev_aln )
                     end
                  end # end score vs maxscore
               end # end isnull
            end # end isvalid
            ridx += 1
            fidx += 1
            continue
         end # is reverse in same gene as fwd
      end # end dist vs. seed_pair_range
      if fwd_sorted[fidx] > rev_sorted[ridx]
         ridx += 1
      else
         fidx += 1
      end
   end #end while

   fwd_res,rev_res
end

ext: jl


text:
using LinearMaps
using Magnus

using PyPlot

function particle_in_box(M = 20,L = 1)
    dx = L/(M+1)
    x = linspace(dx,L-dx,M)
    T = SymTridiagonal(2ones(M), -ones(M-1))/dx^2
    ei = eigfact(T)
    T,ei,x
end

function state(ei, c::Vector, t::Real = 0)
    c /= sqrt(sum(abs2,c))
    sel = 1:length(c)
    c = c.*exp.(-im*ei[:values][sel]*t)
    psi = ei[:vectors][:,sel]*c
    E = dot(ei[:values][sel], abs2.(c))
    psi,E
end

function test_propagation(name, psi0, tmax, period, N, propagator)
    @time psi = integrate(psi0,
                          tmax*period, N,
                          propagator;
                          save_intermediate = true,
                          verbose = true)[:V]
    t = linspace(0,tmax,N+1)
    n = Vector{Float64}(length(t))
    p = Vector{Complex128}(length(t))
    for i = eachindex(t)
        n[i] = norm(psi[:,i])
        p[i] = dot(psi0,psi[:,i])
        psi[:,i] -= p[i]*psi0
    end

    figure(name,figsize=(8,10))
    clf()
    subplot(311)
    pcolormesh(t, x, abs2.(psi), rasterized=true, cmap = plt[:cm][:get_cmap]("viridis"))
    margins(0,0)
    gca()[:set_xticklabels]([])
    title(L"|\Psi-\Psi_0\langle\Psi_0|\Psi\rangle|^2")
    ylabel(L"x")
    subplot(312)
    plot(t,field(t))
    margins(0,0.1)
    gca()[:set_xticklabels]([])
    ylabel(L"F(t)")
    subplot(313)
    plot(t,n, "-")
    margins(0,0.1)
    xlabel(L"t/T")
    ylabel(L"|\Psi|^2")
    a2 = gca()[:twinx]()
    a2[:plot](t,abs.(p), "--")
    margins(0,0.1)
    ylabel(L"|\langle\Psi_0|\Psi\rangle|")
    tight_layout()
end

M = 50
T,ei,x = particle_in_box(M)
psi0,E0 = state(ei, [1])

H₀ = LinearMap(sparse(T))
F₀ = 1e4
period = 1.0/abs(ei[:values][1]-ei[:values][2])
ω = 2π/4
mx = Diagonal(collect(x-mean(x)))
D = LinearMap(sparse(mx))
tmax = 50
field = t -> F₀*exp.(-(t-tmax/2).^2/2dt^2).*sin.(ω*t)
F = t -> field(t*period)
Hᵢ = t -> F(t)*D
H = t -> H₀ + Hᵢ(t)

midpoint = MidpointPropagator(H, -im, LanczosExponentiator(30,psi0))
cfet = CFET4BfCPropagator(H₀, F, D, -im, LanczosExponentiator(30,psi0))

ndt = 300
dt = 6

test_propagation("midpoint", psi0, tmax, period, tmax*300, midpoint)
test_propagation("CFET", psi0, tmax, period, tmax*30, cfet)


text:
function lMap(x)
    μ =  3.9999999999
    return μ*x*(1 - x)
end


text:
"""
An RLEDataFrame extends DataFrame and contains a colection of like-length and like-type
    RLEVectors. In a way, this creates a type like an RLE matrix. But, we deliberately
    avoid the complexity of matrix operations, such as factorization. It is expected
    that most operations will be column-wise. Based on RleDataFrame from Bioconductor's
    `genoset` package (also by Peter Haverty).

### Constructors

```julia
DataFrame(columns::Vector{RLEVector},  names::Vector{Symbol})
DataFrame(kwargs...)
```

### Examples
```julia
RLEDataFrame( [RLEVector([1, 1, 2]),  RLEVector([2, 2, 2])], [:a, :b] )
```

"""
type RLEDataFrame{T1, T2} <: AbstractDataFrame
    columns::RLEVectorList{T1, T2}
    colindex::Index
    function RLEDataFrame(colindex, rvl)
        lens = map(length, rvl)
        (min, max) = extrema(lens)
        if min != max # Redundant with DataFrame
            throw(ArgumentError("All incoming columns must be of equal length."))
        end
        new(rvl, DataFrames.Index(colindex))
    end
end

function RLEDataFrame{T1, T2}(names::Vector{Symbol}, rles::Array{RLEVector{T1, T2}}...)
    rvl = RLEVectorList(rles)
    RLEDataFrame{T1, T2}(names, rvl)
end

function rowSums(df::RLEDataFrame)
    sums = x[1]
    for i in 2:ncol(x)
        sums = sums + x[i]
    end
    return(sums)
end
rowMeans(df::RLEDataFrame) = rowSum(df) ./ ncol(df)
colSums(df::RLEDataFrame) = map(sum, df)
colMeans(df::RLEDataFrame) = colSums(df) ./ nrow(df)


text:
using Makie
using Random
#ϵ ball bridge in likely points with Euler discretization
include("./examples/sin.jl")
#drift(X::LogGrowthSDE, )
drift(X::SinSDE, x) = X.α*sin(x)

function ϵballBridge(S, ϵ, T,discretization, u, v, N)
    Int(T / discretization)
    Ξ = Vector{Float64}[]
    time_grid = 0:discretization:T
    x = fill(0.0, length(time_grid))
    x[1] = u
    while length(Ξ)<N
        for i in 1:length(time_grid)-1
            x[i+1] = x[i] +  drift(S, x[i])*discretization + sqrt(discretization)randn()
        end
        if x[end] < v + ϵ && x[end] > v - ϵ
            push!(Ξ, copy(x))
        end
    end
    return Ξ
end

Random.seed!(1)
ebridge = ϵballBridge(SinSDE(.7,4,5),0.05,50.0,0.0005,-π,π, 500)
#a = plot(0:0.0005:50, ebridge[1:30], leg=false)
#hline!(a, [n*π for n in -3:2:5])
a = Float64[]
push!(a, 1.0)
function scatter_joint(a, ebridge)
    x1 = Float64[]
    x2 = Float64[]
    for i in 1:length(ebridge)
        push!(x1, ebridge[i][Int((length(ebridge[i]) - 1)/4) + 1])
        push!(x2, ebridge[i][Int((length(ebridge[i]) - 1)*3/4) + 1])
    end
    Makie.scatter!(a, x1, x2, color = (:darkred, 1.0), markersize = 0.25)
end
using Plots
y = runall(true)

l1 = Λ(50/4, 0, 0, 50)
l2 = Λ(50*3/4, 1, 1, 50)

y1 = [1/4*π + 3/4*(-π) + y[i].ξ[1]*l1 + y[i].ξ[2]*l2 for i in 1:length(y)]
y2 = [1/4*(-π) + 3/4*π + y[i].ξ[1]*l1 + y[i].ξ[3]*l2 for i in 1:length(y)]

using Makie
ed = lines(y1, y2, color = (:lightblue), linewidth=0.5)
lines!(ed, y1, y2, color = (:blue, 0.01), linewidth=1.0)
scatter_joint(ed, ebridge )
axis = ed[Axis]
axis[:names, :axisnames] = ("X_12.5", "X_37.5")
Makie.save("joint_distribution.png", ed)

Makie.scatter!

max_stars_repo_name: SebaGraz/ZZDiffusionBridge

max_issues_repo_name: SebaGraz/ZZDiffusionBridge

max_forks_repo_name: SebaGraz/ZZDiffusionBridge


text:
@doc raw"""
    ProbabilitySimplex{n} <: AbstractEmbeddedManifold{ℝ,DefaultEmbeddingType}

The (relative interior of) the probability simplex is the set
````math
Δ^n := \biggl\{ p ∈ ℝ^{n+1}\ \big|\ p_i > 0 \text{ for all } i=1,…,n+1,
\text{ and } ⟨\mathbb{1},p⟩ = \sum_{i=1}^{n+1} p_i = 1\biggr\},
````
where $\mathbb{1}=(1,…,1)^{\mathrm{T}}∈ ℝ^{n+1}$ denotes the vector containing only ones.

This set is also called the unit simplex or standard simplex.

The tangent space is given by
````math
T_pΔ^n = \biggl\{ X ∈ ℝ^{n+1}\ \big|\ ⟨\mathbb{1},X⟩ = \sum_{i=1}^{n+1} X_i = 0 \biggr\}
````

The manifold is implemented assuming the Fisher-Rao metric for the multinomial distribution,
which is equivalent to the induced metric from isometrically embedding the probability
simplex in the $n$-sphere of radius 2.
The corresponding diffeomorphism $\varphi: \mathbb Δ^n → \mathcal N$,
where $\mathcal N \subset 2𝕊^n$ is given by $\varphi(p) = 2\sqrt{p}$.

This implementation follows the notation in [^ÅströmPetraSchmitzerSchnörr2017].

[^ÅströmPetraSchmitzerSchnörr2017]:
    > F. Åström, S. Petra, B. Schmitzer, C. Schnörr: “Image Labeling by Assignment”,
    > Journal of Mathematical Imaging and Vision, 58(2), pp. 221–238, 2017.
    > doi: [10.1007/s10851-016-0702-4](https://doi.org/10.1007/s10851-016-0702-4)
    > arxiv: [1603.05285](https://arxiv.org/abs/1603.05285).
"""
struct ProbabilitySimplex{n} <: AbstractEmbeddedManifold{ℝ,DefaultEmbeddingType} end

ProbabilitySimplex(n::Int) = ProbabilitySimplex{n}()

"""
    SoftmaxRetraction <: AbstractRetractionMethod

Describes a retraction that is based on the softmax function.
"""
struct SoftmaxRetraction <: AbstractRetractionMethod end

"""
    SoftmaxInverseRetraction <: AbstractInverseRetractionMethod

Describes an inverse retraction that is based on the softmax function.
"""
struct SoftmaxInverseRetraction <: AbstractInverseRetractionMethod end

"""
    FisherRaoMetric <: Metric

The Fisher-Rao metric or Fisher information metric is a particular Riemannian metric which
can be defined on a smooth statistical manifold, i.e., a smooth manifold whose points are
probability measures defined on a common probability space.

See for example the [`ProbabilitySimplex`](@ref).
"""
struct FisherRaoMetric <: Metric end

"""
    check_manifold_point(M::ProbabilitySimplex, p; kwargs...)

Check whether `p` is a valid point on the [`ProbabilitySimplex`](@ref) `M`, i.e. is a point in
the embedding with positive entries that sum to one
The tolerance for the last test can be set using the `kwargs...`.
"""
function check_manifold_point(M::ProbabilitySimplex, p; kwargs...)
    mpv = invoke(
        check_manifold_point,
        Tuple{(typeof(get_embedding(M))),typeof(p)},
        get_embedding(M),
        p;
        kwargs...,
    )
    mpv === nothing || return mpv
    if minimum(p) <= 0
        return DomainError(
            minimum(p),
            "The point $(p) does not lie on the $(M) since it has nonpositive entries.",
        )
    end
    if !isapprox(sum(p), 1.0; kwargs...)
        return DomainError(
            sum(p),
            "The point $(p) does not lie on the $(M) since its sum is not 1.",
        )
    end
    return nothing
end

"""
    check_tangent_vector(M::ProbabilitySimplex, p, X; check_base_point = true, kwargs... )

Check whether `X` is a tangent vector to `p` on the [`ProbabilitySimplex`](@ref) `M`, i.e.
after [`check_manifold_point`](@ref check_manifold_point(::ProbabilitySimplex, ::Any))`(M,p)`,
`X` has to be of same dimension as `p` and its elements have to sum to one.
The optional parameter `check_base_point` indicates, whether to call
[`check_manifold_point`](@ref check_manifold_point(::ProbabilitySimplex, ::Any))  for `p` or not.
The tolerance for the last test can be set using the `kwargs...`.
"""
function check_tangent_vector(M::ProbabilitySimplex, p, X; check_base_point=true, kwargs...)
    if check_base_point
        mpe = check_manifold_point(M, p; kwargs...)
        mpe === nothing || return mpe
    end
    mpv = invoke(
        check_tangent_vector,
        Tuple{typeof(get_embedding(M)),typeof(p),typeof(X)},
        get_embedding(M),
        p,
        X;
        check_base_point=false, # already checked above
        kwargs...,
    )
    mpv === nothing || return mpv
    if !isapprox(sum(X), 0.0; kwargs...)
        return DomainError(
            sum(X),
            "The vector $(X) is not a tangent vector to $(p) on $(M), since its elements to not sum up to 0.",
        )
    end
    return nothing
end

decorated_manifold(M::ProbabilitySimplex) = Euclidean(representation_size(M)...; field=ℝ)

default_metric_dispatch(::ProbabilitySimplex, ::FisherRaoMetric) = Val(true)

@doc raw"""
    distance(M,p,q)

Compute the distance between two points on the [`ProbabilitySimplex`](@ref) `M`.
The formula reads
````math
d_{Δ^n}(p,q) = 2\arccos \biggl( \sum_{i=1}^{n+1} \sqrt{p_i q_i} \biggr)
````
"""
function distance(::ProbabilitySimplex, p, q)
    sumsqrt = zero(Base.promote_eltype(p, q))
    @inbounds for i in eachindex(p, q)
        sumsqrt += sqrt(p[i] * q[i])
    end
    return 2 * acos(sumsqrt)
end

@doc raw"""
    exp(M::ProbabilitySimplex,p,X)

Compute the exponential map on the probability simplex.

````math
\exp_pX = \frac{1}{2}\Bigl(p+\frac{X_p^2}{\lVert X_p \rVert^2}\Bigr)
+ \frac{1}{2}\Bigl(p - \frac{X_p^2}{\lVert X_p \rVert^2}\Bigr)\cos(\lVert X_p\rVert)
+ \frac{1}{\lVert X_p \rVert}\sqrt{p}\sin(\lVert X_p\rVert),
````

where $X_p = \frac{X}{\sqrt{p}}$, with its division meant elementwise, as well as for the
operations $X_p^2$ and $\sqrt{p}$.
"""
exp(::ProbabilitySimplex, ::Any...)

function exp!(::ProbabilitySimplex, q, p, X)
    s = sqrt.(p)
    Xs = X ./ s ./ 2
    θ = norm(Xs)
    q .= (cos(θ) .* s .+ usinc(θ) .* Xs) .^ 2
    return q
end

@doc raw"""
    injectivity_radius(M,p)

compute the injectivity radius on the [`ProbabilitySimplex`](@ref) `M` at the point `p`,
i.e. the distanceradius to a point near/on the boundary, that could be reached by following the
geodesic.
"""
function injectivity_radius(::ProbabilitySimplex{n}, p) where {n}
    i = argmin(p)
    s = sum(p) - p[i]
    return 2 * acos(sqrt(s))
end
function injectivity_radius(M::ProbabilitySimplex, p, ::ExponentialRetraction)
    return injectivity_radius(M, p)
end
injectivity_radius(M::ProbabilitySimplex, p, ::SoftmaxRetraction) = injectivity_radius(M, p)
injectivity_radius(M::ProbabilitySimplex) = 0
injectivity_radius(M::ProbabilitySimplex, ::SoftmaxRetraction) = 0
injectivity_radius(M::ProbabilitySimplex, ::ExponentialRetraction) = 0
eval(
    quote
        @invoke_maker 1 Manifold injectivity_radius(
            M::ProbabilitySimplex,
            rm::AbstractRetractionMethod,
        )
    end,
)

@doc raw"""
    inner(M::ProbabilitySimplex,p,X,Y)

Compute the inner product of two tangent vectors `X`, `Y` from the tangent space $T_pΔ^n$ at
`p`. The formula reads
````math
g_p(X,Y) = \sum_{i=1}^{n+1}\frac{X_iY_i}{p_i}
````
"""
function inner(::ProbabilitySimplex, p, X, Y)
    d = zero(Base.promote_eltype(p, X, Y))
    @inbounds for i in eachindex(p, X, Y)
        d += X[i] * Y[i] / p[i]
    end
    return d
end

@doc raw"""
    inverse_retract(M::ProbabilitySimplex, p, q, ::SoftmaxInverseRetraction)

Compute a first order approximation by projection. The formula reads
````math
\operatorname{retr}^{-1}_p q = \bigl( I_{n+1} - \frac{1}{n}\mathbb{1}^{n+1,n+1} \bigr)(\log(q)-\log(p))
````
where $\mathbb{1}^{m,n}$ is the size `(m,n)` matrix containing ones, and $\log$ is applied elementwise.
"""
inverse_retract(::ProbabilitySimplex, ::Any, ::Any, ::SoftmaxInverseRetraction)

function inverse_retract!(
    ::ProbabilitySimplex{n},
    X,
    p,
    q,
    ::SoftmaxInverseRetraction,
) where {n}
    X .= log.(q) .- log.(p)
    meanlogdiff = mean(X)
    X .-= meanlogdiff
    return X
end

@doc raw"""
    log(M::ProbabilitySimplex, p, q)

Compute the logarithmic map of `p` and `q` on the [`ProbabilitySimplex`](@ref) `M`.

````math
\log_pq = \frac{d_{Δ^n}(p,q)}{\sqrt{1-⟨\sqrt{p},\sqrt{q}⟩}}(\sqrt{pq} - ⟨\sqrt{p},\sqrt{q}⟩p),
````

where $pq$ and $\sqrt{p}$ is meant elementwise.
"""
log(::ProbabilitySimplex, ::Any...)

function log!(::ProbabilitySimplex, X, p, q)
    if p ≈ q
        fill!(X, 0)
    else
        z = sqrt.(p .* q)
        s = sum(z)
        X .= 2 * acos(s) / sqrt(1 - s^2) .* (z .- s .* p)
    end
    return X
end

@doc raw"""
    manifold_dimension(M::ProbabilitySimplex{n})

Returns the manifold dimension of the probability simplex in $ℝ^{n+1}$, i.e.
````math
    \dim_{Δ^n} = n.
````
"""
manifold_dimension(::ProbabilitySimplex{n}) where {n} = n

@doc raw"""
    mean(
        M::ProbabilitySimplex,
        x::AbstractVector,
        [w::AbstractWeights,]
        method = GeodesicInterpolation();
        kwargs...,
    )

Compute the Riemannian [`mean`](@ref mean(M::Manifold, args...)) of `x` using
[`GeodesicInterpolation`](@ref).
"""
mean(::ProbabilitySimplex, ::Any...)

function Statistics.mean!(
    M::ProbabilitySimplex,
    p,
    x::AbstractVector,
    w::AbstractVector;
    kwargs...,
)
    return mean!(M, p, x, w, GeodesicInterpolation(); kwargs...)
end

@doc raw"""
    project(M::ProbabilitySimplex, p, Y)

project `Y` from the embedding onto the tangent space at `p` on
the [`ProbabilitySimplex`](@ref) `M`. The formula reads

````math
\operatorname{proj}_{Δ^n}(p,Y) = Y - ⟨p,Y⟩p.
````
"""
project(::ProbabilitySimplex, ::Any, ::Any)

function project!(::ProbabilitySimplex, q, p)
    if any(x -> x <= 0, p)
        throw(
            DomainError(
                p,
                "All coordinates of point from the embedding, that should be projected, must be positive, otherwise the projection is not well defined.",
            ),
        )
    end
    q .= p ./ sum(p)
    return q
end

function project!(::ProbabilitySimplex, X, p, Y)
    X .= Y .- sum(Y) .* p
    return X
end

@doc raw"""
    representation_size(::ProbabilitySimplex{n})

return the representation size of points in the $n$-dimensional probability simplex,
i.e. an array size of `(n+1,)`.
"""
representation_size(::ProbabilitySimplex{n}) where {n} = (n + 1,)

@doc raw"""
    retract(M::ProbabilitySimplex, p, X, ::SoftmaxRetraction)

Compute a first order approximation by applying the softmax function. The formula reads

````math
\operatorname{retr}_p X = \frac{p\mathrm{e}^X}{⟨p,\mathrm{e}^X⟩},
````

where multiplication, exponentiation and division are meant elementwise.
"""
retract(::ProbabilitySimplex, ::Any, ::Any, ::SoftmaxRetraction)

function retract!(::ProbabilitySimplex, q, p, X, ::SoftmaxRetraction)
    s = zero(eltype(q))
    @inbounds for i in eachindex(q, p, X)
        q[i] = p[i] * exp(X[i])
        s += q[i]
    end
    q ./= s
    return q
end

function Base.show(io::IO, ::ProbabilitySimplex{n}) where {n}
    return print(io, "ProbabilitySimplex($(n))")
end

@doc raw"""
    zero_tangent_vector(M::ProbabilitySimplex,p)

returns the zero tangent vector in the tangent space of the point `p`  from the
[`ProbabilitySimplex`](@ref) `M`, i.e. its representation by the zero vector in the embedding.
"""
zero_tangent_vector(::ProbabilitySimplex, ::Any)

zero_tangent_vector!(M::ProbabilitySimplex, v, p) = fill!(v, 0)


text:
using MolecularGraph:
    fgroupquery, fgroupcond, fgrouprecord, getterm, setterm!

@testset "funcgroup" begin

@testset "fgroupquery" begin
    glyc = smilestomol("C(O)C(O)CO")
    res = fgroupquery(glyc, "[#6][OD1]")
    @test length(res) == 3
end

@testset "fgroupcond" begin
    alkoxide = smilestomol("C(O)C(O)C[O-]")
    res = fgroupcond(alkoxide, Dict("any" => ["[#6][OD1]", "[#6][O-]"]))
    @test length(res) == 3
end

@testset "fgrouprecord" begin
    coumarin = smilestomol("C1=CC(=O)OC2=C1C=CC=C2")
    fgc = FunctionalGroupClassifier()
    setterm!(fgc, :Carbonyl, fgroupquery(coumarin, "[#6]=[OD1]"))
    res1 = fgrouprecord(
        coumarin, fgc, Dict("have" => ["Carbonyl"], "query" => "[#6](=O)[#8]"))
    @test issetequal(collect(res1)[1], [3, 4, 5])

    res2 = fgrouprecord(
        coumarin, fgc,
        Dict("have" => ["Carbonyl", "Something"], "query" => "[#6](=O)[#8]"))
    @test isempty(res2)

    setterm!(fgc, :Carboxyl, res1)
    res3 = fgrouprecord(
        coumarin, fgc, Dict("isa" => ["Carboxyl"], "query" => "c(=O)o"))
    @test issetequal(collect(res3)[1], [3, 4, 5])
    res4 = fgrouprecord(
        coumarin, fgc, Dict("isa" => ["Carboxyl"], "query" => "[#6](=O)O"))
    @test isempty(res4)
end

@testset "functionalgroupgraph" begin
    paclitaxel = smilestomol(
        "CC1=C2[C@@]([C@]([C@H]([C@@H]3[C@]4([C@H](OC4)C[C@@H]([C@]3(C(=O)[C@@H]2OC(=O)C)C)O)OC(=O)C)OC(=O)c5ccccc5)(C[C@@H]1OC(=O)[C@H](O)[C@@H](NC(=O)c6ccccc6)c7ccccc7)O)(C)C"
    )
    paclitaxel = removehydrogens(paclitaxel)
    fgc = functionalgroupgraph(paclitaxel)
    @test length(getterm(fgc, :Carboxyl)) == 4
    @test length(getterm(fgc, :SecAlcohol)) == 2
    @test length(getterm(fgc, :FourMemberedRing)) == 1
    @test length(getterm(fgc, :EightMemberedRing)) == 1
end

end # funcgroup


text:
mutable struct FedProxServer{T1<:Int64, T2<:Float64, T3<:SparseMatrixCSC{Float64, Int64}, T4<:Vector{Int64}, T5<:Vector{FedProxClient}, T6<:Matrix{Float64}} <:AbstractServer
    Xtest::T3                        # testing data
    Ytest::T4                        # testing label
    num_classes::T1                  # number of classes
    num_clients::T1                  # number of clients
    participation_rate::T2           # participation rate
    clients::T5                      # set of clients
    W::T6                           # global models
    τ::T1                           # number of particiating clients
    selectedIndices::T4             # selected clients
    function FedProxServer(Xtest::SparseMatrixCSC{Float64, Int64}, Ytest::Vector{Int64}, config::Dict{String, Real})
        num_classes = config["num_classes"]
        num_clients = config["num_clients"]
        participation_rate = config["participation_rate"]
        τ = floor(Int64, num_clients * participation_rate)
        _, d = size(Xtest)
        clients = Vector{FedProxClient}(undef, num_clients)
        W = zeros(Float64, d, num_classes)
        new{Int64, Float64, SparseMatrixCSC{Float64, Int64}, Vector{Int64}, Vector{FedProxClient}, Matrix{Float64}}(Xtest, Ytest, num_classes, num_clients, participation_rate, clients, W, τ)
    end
end

# Select clients
function select!(
    server::FedProxServer
    )
    server.selectedIndices = randperm(server.num_clients)[1:server.τ]
    return nothing
end

# Send global model
function sendModel!(
    server::FedProxServer
)
    # Only send model to selected clients
    for idx in server.selectedIndices
        server.clients[idx].W = copy(server.W)
    end
    return nothing
end

# Send global model to all clients
function sendModelToAllClients!(
    server::FedProxServer
)
    # Only send model to selected clients
    for client in server.clients
        client.W = copy(server.W)
    end
    return nothing
end

# Aggregation
function aggregate!(
    server::FedProxServer
)
    # Take the average of the selected clients' model
    fill!(server.W, 0.0)
    for i = 1:server.τ
        idx = server.selectedIndices[i]
        server.W += server.clients[idx].W
    end
    server.W ./= server.τ
    return nothing
end

# Calculate objective value
function getObjValue(
    server::FedProxServer
)
    objValue = 0.0
    totalNumData = 0
    for i = 1:server.num_clients
        n = size(server.clients[i].Xtrain, 1)
        totalNumData += n
        objValue += n * getObjValue(server.clients[i])
    end
    return objValue / totalNumData
end

########################################################################################################
mutable struct FedProxServerNN{T1<:Int64, T2<:Float64, T3<:SparseMatrixCSC{Float64, Int64}, T4<:Flux.OneHotArray, T5<:Vector{FedProxClientNN}, T6<:Flux.Chain, T7<:Vector{Int64}} <:AbstractServer
    Xtest::T3                        # testing data
    Ytest::T4                        # testing label
    num_classes::T1                  # number of classes
    num_clients::T1                  # number of clients
    participation_rate::T2           # participation rate
    clients::T5                      # set of clients
    W::T6                           # global models
    τ::T1                           # number of particiating clients
    selectedIndices::T7             # selected clients
    function FedProxServerNN(Xtest::SparseMatrixCSC{Float64, Int64}, Ytest::Vector{Int64}, dim::Int64, config::Dict{String, Real})
        num_classes = config["num_classes"]
        Ytest = Flux.onehotbatch(Ytest, 1:num_classes)
        num_clients = config["num_clients"]
        participation_rate = config["participation_rate"]
        W = Chain( Dense(780, dim, relu, bias=false), Dense(dim, 10, bias=false), NNlib.softmax)
        τ = floor(Int64, num_clients * participation_rate)
        clients = Vector{FedProxClientNN}(undef, num_clients)
        new{Int64, Float64, SparseMatrixCSC{Float64, Int64}, Flux.OneHotArray, Vector{FedProxClientNN}, Flux.Chain, Vector{Int64}}(Xtest, Ytest, num_classes, num_clients, participation_rate, clients, W, τ)
    end
end

# Select clients
function select!(
    server::FedProxServerNN
    )
    server.selectedIndices = randperm(server.num_clients)[1:server.τ]
    return nothing
end

# Send global model, now exactly the same as `synchronize!`.
function sendModel!(
    server::FedProxServerNN
)
    # Only send model to selected clients
    for i in server.selectedIndices
        c = server.clients[i]
        for j = 1:length(params(server.W))
            params(c.W)[j] .= params(server.W)[j]
        end
    end
    return nothing
end

# Send global model to all clients
function sendModelToAllClients!(
    server::FedProxServerNN
)
    # Only send model to selected clients
    for c in server.clients
        for j = 1:length(params(server.W))
            params(c.W)[j] .= params(server.W)[j]
        end
    end
    return nothing
end

# Aggregation
function aggregate!(
    server::FedProxServerNN
)
    # Take the average of the selected clients' model
    l = length(params(server.W))
    for j = 1:l
        fill!(params(server.W)[j], 0.0)
    end
    for i = 1:server.τ
        idx = server.selectedIndices[i]
        c = server.clients[idx]
        for j = 1:l
            params(server.W)[j] .+= (1/server.τ) * params(c.W)[j]
        end
    end
    return nothing
end

# synchronization
function synchronize!(
    server::FedProxServerNN
)
    for i in server.selectedIndices
        c = server.clients[i]
        for j = 1:length(params(server.W))
            params(c.W)[j] .= params(server.W)[j]
        end
    end
    return nothing
end

# Calculate objective value
function getObjValue(
    server::FedProxServerNN
)
    objValue = 0.0
    totalNumData = 0
    for i = 1:server.num_clients
        c = server.clients[i]
        n = size(c.Xtrain, 1)
        totalNumData += n
        objValue += n * obj( c.XtrainT, c.Ytrain, server.W, c.λ)
    end
    return objValue / totalNumData
end