# using LatticeGPU
using ADerrors
using BDIO


function read_ff(name::String)

    file = BDIO_open(name,"r","FerFlow correlators")

    while BDIO_get_uinfo(file) != 14
        BDIO_seek!(file)
    end
    ihdr = Vector{Int32}(undef,1)
    BDIO_read(file, ihdr)
    if ihdr[1] == 1730280201
        println("Reading file ",name)
    else
        error("Wrong IHDR in file ", name)
    end

    while BDIO_get_uinfo(file) != 1
        BDIO_seek!(file)
    end

    dims=Vector{Int32}(undef,1)
    iL=Vector{Int32}(undef,4)
    N_noiseff=Vector{Int32}(undef,1)
    nsteps=Vector{Int32}(undef,1)
    N_noisebf=Vector{Int32}(undef,1)
    nflow=Vector{Int32}(undef,1)
    tzeroff=Vector{Float64}(undef,1)
    epsilon=Vector{Float64}(undef,1)
    tsourceff=Vector{Int32}(undef,1)
    tsourcebf=Vector{Int32}(undef,1)
    Nf=Vector{Int32}(undef,1)


    BDIO_read(file, dims)
    BDIO_read(file, iL)
    BDIO_read(file, Nf)

    dpars=[Vector{Float64}(undef,2) for _ in 1:Nf[1]] # m0,csw
    th=[Vector{ComplexF64}(undef,4) for _ in 1:Nf[1]]

    for i in 1:Nf[1]
        BDIO_read(file, dpars[i])
        BDIO_read(file, th[i])
    end

    BDIO_read(file, N_noiseff)
    BDIO_read(file, nsteps)
    BDIO_read(file, N_noisebf)
    BDIO_read(file, nflow)
    BDIO_read(file, tzeroff)
    BDIO_read(file, epsilon)
    BDIO_read(file, tsourceff)
    BDIO_read(file, tsourcebf)
    flow_times=Vector{Float64}(undef,nflow[1])
    nflow[1] > 0 ? BDIO_read(file, flow_times) : nothing

    println("Reading data from a ",Int64.(iL), " lattice.")
    println("Frontflow: ",N_noiseff[1], " sources at euclidean time ", tsourceff[1] ,". Flowtime starting at ",tzeroff[1], " with ",nsteps[1]," steps of size ", epsilon[1])
    println("Backflow: ",N_noisebf[1]," sources at euclidean time ",tsourcebf[1],". Flow times are ", flow_times)

    N_noiseff = N_noiseff[1]
    N_noisebf = N_noisebf[1]
    nflow = nflow[1]
    nsteps = nsteps[1]
    Nf = Nf[1]

    TvecR = Vector{Float64}(undef,iL[4]);
    TvecC = Vector{ComplexF64}(undef,iL[4]);
    CN = Vector{ComplexF64}(undef,1);
    RN = Vector{Float64}(undef,1);

    pp_corr = Array{Float64}(undef,iL[4],N_noiseff,Nf);
    ap_corr = Array{ComplexF64}(undef,iL[4],N_noiseff,Nf);
    pphat0 = Array{Float64}(undef,iL[4],N_noiseff,Nf);
    pptilde0 = Array{ComplexF64}(undef,iL[4],N_noiseff,Nf);
    pp_corr_t = Array{Float64}(undef,iL[4],N_noiseff,nsteps+1,Nf);
    ap_corr_t = Array{ComplexF64}(undef,iL[4],N_noiseff,nsteps+1,Nf);
    pphat_t = Array{Float64}(undef,iL[4],N_noiseff,nsteps+1,Nf);
    pptilde_t = Array{ComplexF64}(undef,iL[4],N_noiseff,nsteps+1,Nf);
    Eoft = Array{Complex{Float64}}(undef,1+nsteps,iL[4]);
    Eoft0 = Array{Complex{Float64}}(undef,iL[4]);
    Qt = Array{Complex{Float64}}(undef,1+nsteps,iL[4]);
    Qt0 = Array{Complex{Float64}}(undef,iL[4]);

    pp_corr_tfl = Array{Float64}(undef,iL[4],N_noisebf,nflow,Nf);
    ap_corr_tfl = Array{ComplexF64}(undef,iL[4],N_noisebf,nflow,Nf);
    Sigma = Array{ComplexF64}(undef,iL[4],N_noisebf,nflow,Nf);
    Sigma_cfl = Array{ComplexF64}(undef,iL[4],N_noisebf,nflow,Nf);
    Sigma2 = Array{ComplexF64}(undef,iL[4],N_noisebf,nflow,Nf);
    Sigma2_cfl = Array{ComplexF64}(undef,iL[4],N_noisebf,nflow,Nf);
    ChiDchi = Array{ComplexF64}(undef,iL[4],N_noisebf,nflow,Nf);
    ChiDchi_cfl = Array{ComplexF64}(undef,iL[4],N_noisebf,nflow,Nf);

    pp_corr_mc = [Vector{Float64}() for _ in 1:iL[4], _ in 1:N_noiseff, _ in 1:Nf]
    ap_corr_mc = [Vector{ComplexF64}() for _ in 1:iL[4], _ in 1:N_noiseff, _ in 1:Nf]
    pphat0_mc = [Vector{Float64}() for _ in 1:iL[4], _ in 1:N_noiseff, _ in 1:Nf]
    pptilde0_mc = [Vector{ComplexF64}() for _ in 1:iL[4], _ in 1:N_noiseff, _ in 1:Nf]
    pp_corr_t_mc = [Vector{Float64}() for _ in 1:iL[4], _ in 1:N_noiseff, _ in 1:nsteps+1, _ in 1:Nf]
    ap_corr_t_mc = [Vector{ComplexF64}() for _ in 1:iL[4], _ in 1:N_noiseff, _ in 1:nsteps+1, _ in 1:Nf]
    pphat_t_mc = [Vector{Float64}() for _ in 1:iL[4], _ in 1:N_noiseff, _ in 1:nsteps+1, _ in 1:Nf]
    pptilde_t_mc = [Vector{ComplexF64}() for _ in 1:iL[4], _ in 1:N_noiseff, _ in 1:nsteps+1, _ in 1:Nf]
    Eoft_mc = [Vector{Complex{Float64}}() for _ in 1:1+nsteps, _ in 1:iL[4]]
    Eoft0_mc = [Vector{Complex{Float64}}() for _ in 1:iL[4]]
    Qt_mc = [Vector{Complex{Float64}}() for _ in 1:1+nsteps, _ in 1:iL[4]]
    Qt0_mc = [Vector{Complex{Float64}}() for _ in 1:iL[4]]

    pp_corr_tfl_mc = [Vector{Float64}() for _ in 1:iL[4], _ in 1:N_noisebf, _ in 1:nflow, _ in 1:Nf]
    ap_corr_tfl_mc = [Vector{ComplexF64}() for _ in 1:iL[4], _ in 1:N_noisebf, _ in 1:nflow, _ in 1:Nf]
    Sigma_mc = [Vector{ComplexF64}() for _ in 1:iL[4], _ in 1:N_noisebf, _ in 1:nflow, _ in 1:Nf]
    Sigma_cfl_mc = [Vector{ComplexF64}() for _ in 1:iL[4], _ in 1:N_noisebf, _ in 1:nflow, _ in 1:Nf]
    Sigma2_mc = [Vector{ComplexF64}() for _ in 1:iL[4], _ in 1:N_noisebf, _ in 1:nflow, _ in 1:Nf]
    Sigma2_cfl_mc = [Vector{ComplexF64}() for _ in 1:iL[4], _ in 1:N_noisebf, _ in 1:nflow, _ in 1:Nf]
    ChiDchi_mc = [Vector{ComplexF64}() for _ in 1:iL[4], _ in 1:N_noisebf, _ in 1:nflow, _ in 1:Nf]
    ChiDchi_cfl_mc = [Vector{ComplexF64}() for _ in 1:iL[4], _ in 1:N_noisebf, _ in 1:nflow, _ in 1:Nf]

    while BDIO_get_uinfo(file) != 8
        BDIO_seek!(file)
    end

    while BDIO_get_uinfo(file) == 8

        for f in 1:Nf
            for noi in 1:N_noiseff
                BDIO_read(file,TvecR)
                pp_corr[:,noi,f] .= TvecR

                BDIO_read(file,TvecC)
                ap_corr[:,noi,f] .= TvecC

                BDIO_read(file,TvecR)
                pphat0[:,noi,f] .= TvecR

                BDIO_read(file,TvecC)
                pptilde0[:,noi,f] .= TvecC

                for fl in 1:nsteps+1
                    BDIO_read(file,TvecR)
                    pp_corr_t[:,noi,fl,f].= TvecR

                    BDIO_read(file,TvecC)
                    ap_corr_t[:,noi,fl,f].= TvecC

                    BDIO_read(file,TvecR)
                    pphat_t[:,noi,fl,f].= TvecR

                    BDIO_read(file,TvecC)
                    pptilde_t[:,noi,fl,f].= TvecC
                end
            end
        end

        BDIO_read(file,TvecC)
        Eoft0[:] .= TvecC
        BDIO_read(file,TvecC)
        Qt0[:] .= TvecC

        for fl in 1:nsteps+1
            BDIO_read(file,TvecC)
            Eoft[fl,:] .= TvecC
            BDIO_read(file,TvecC)
            Qt[fl,:] .= TvecC
        end

        for f in 1:Nf
            for noi in 1:N_noisebf
                for fl in 1:nflow
                    BDIO_read(file,TvecR)
                    pp_corr_tfl[:,noi,fl,f].= TvecR

                    BDIO_read(file,TvecC)
                    ap_corr_tfl[:,noi,fl,f].= TvecC

                    BDIO_read(file,TvecC)
                    Sigma[:,noi,fl,f] .= TvecC

                    BDIO_read(file,TvecC)
                    Sigma_cfl[:,noi,fl,f] .= TvecC

                    BDIO_read(file,TvecC)
                    Sigma2[:,noi,fl,f] .= TvecC

                    BDIO_read(file,TvecC)
                    Sigma2_cfl[:,noi,fl,f] .= TvecC

                    BDIO_read(file,TvecC)
                    ChiDchi[:,noi,fl,f] .= TvecC

                    BDIO_read(file,TvecC)
                    ChiDchi_cfl[:,noi,fl,f] .= TvecC
                end
            end
        end

        push!.(Eoft0_mc, Eoft0)
        push!.(Eoft_mc, Eoft)
        push!.(Qt0_mc, Qt0)
        push!.(Qt_mc, Qt)

        push!.(pp_corr_mc, pp_corr)
        push!.(ap_corr_mc, ap_corr)
        push!.(pphat0_mc, pphat0)
        push!.(pptilde0_mc, pptilde0)
        push!.(pp_corr_t_mc, pp_corr_t)
        push!.(ap_corr_t_mc, ap_corr_t)
        push!.(pphat_t_mc, pphat_t)
        push!.(pptilde_t_mc, pptilde_t)

        push!.(pp_corr_tfl_mc, pp_corr_tfl)
        push!.(ap_corr_tfl_mc, ap_corr_tfl)
        push!.(Sigma_mc, Sigma)
        push!.(Sigma_cfl_mc, Sigma_cfl)
        push!.(Sigma2_mc, Sigma2)
        push!.(Sigma2_cfl_mc, Sigma2_cfl)
        push!.(ChiDchi_mc, ChiDchi)
        push!.(ChiDchi_cfl_mc, ChiDchi_cfl)

        BDIO_seek!(file)
        BDIO_seek!(file)
    end

    BDIO_close!(file)

    return Eoft0_mc,Eoft_mc,Qt0_mc,Qt_mc,pp_corr_mc,ap_corr_mc,pp_corr_t_mc,ap_corr_t_mc,pp_corr_tfl_mc,ap_corr_tfl_mc,Sigma_mc,Sigma_cfl_mc,Sigma2_mc,Sigma2_cfl_mc,ChiDchi_mc,ChiDchi_cfl_mc,pphat_t_mc,pptilde_t_mc,pphat0_mc,pptilde0_mc
end


"""
function uwff(file::String)

Loads in the variables 'Eoft','Eoft0','Qt','Qt0', 'pp_corr', 'ap_corr', 'pp_corr_t','ap_corr_t', 'pphat_t', 'pptilde_t','pphat0', 'pptilde0','pp_corr_tfl','ap_corr_tfl' 'Sigma'(v1 and v2), 'Sigma_cfl'(v1 and v2), 'ChiDchi' and 'ChiDchi_cfl' the corresponding
'uwreal' quantities. Noise average has been performed. The indices are, when corresponding, eucliden time and flow times
"""
function uwff(file::String)
    Eoft0_mc,Eoft_mc,Qt0_mc,Qt_mc,pp_corr_mc,ap_corr_mc,pp_corr_t_mc,ap_corr_t_mc,pp_corr_tfl_mc,ap_corr_tfl_mc,Sigma_mc,Sigma_cfl_mc,Sigma2_mc,Sigma2_cfl_mc,ChiDchi_mc,ChiDchi_cfl_mc,pphat_t_mc,pptilde_t_mc,pphat0_mc,pptilde0_mc = read_ff(file)
    runame = String(split(split(file,"/")[end],".")[1])

    T = size(pp_corr_mc)[1]
    Nsff = size(pp_corr_mc)[2]
    nsteps = size(pp_corr_t_mc)[3]-1
    Nsbf = size(pp_corr_tfl_mc)[2]
    Nfl = size(pp_corr_tfl_mc)[3]
    Nf = size(pp_corr_mc)[3]

    global Eoft = [uwreal(real.(Eoft_mc[i,j]) ,runame) for j in 1:T, i in 1:nsteps+1]
    global Eoft0 = [uwreal(real.(Eoft0_mc[j]) ,runame) for j in 1:T]

    global Qt = [uwreal(real.(Qt_mc[i,j]) ,runame) for j in 1:T, i in 1:nsteps+1]
    global Qt0 = [uwreal(real.(Qt0_mc[j]) ,runame) for j in 1:T]

    pp_corr_mc = sum(pp_corr_mc,dims = 2)[:,1,:]./Nsff
    global pp_corr = [uwreal(pp_corr_mc[i,f] ,runame) for i in 1:T, f in 1:Nf]

    ap_corr_mc = real.(sum(ap_corr_mc,dims = 2)[:,1,:]./Nsff)
    global ap_corr = [uwreal(ap_corr_mc[i,f] ,runame) for i in 1:T, f in 1:Nf]

    global pphat0 = [uwreal(sum(pphat0_mc,dims = 2)[i,1,f]./Nsff ,runame) for i in 1:T, f in 1:Nf]
    global pptilde0 = [uwreal(real.(sum(pptilde0_mc,dims = 2)[i,1,f]./Nsff) ,runame) for i in 1:T, f in 1:Nf]

    pp_corr_t_mc = sum(pp_corr_t_mc,dims = 2)[:,1,:,:]./Nsff
    global pp_corr_t = [uwreal(pp_corr_t_mc[i,k,f] ,runame) for i in 1:T, k in 1:nsteps+1, f in 1:Nf]
    ap_corr_t_mc = real.(sum(ap_corr_t_mc,dims = 2)[:,1,:,:]./Nsff)
    global ap_corr_t = [uwreal(ap_corr_t_mc[i,k,f] ,runame) for i in 1:T, k in 1:nsteps+1, f in 1:Nf]

    pphat_t_mc = sum(pphat_t_mc,dims = 2)[:,1,:,:]./Nsff
    global pphat_t = [uwreal(pphat_t_mc[i,k,f] ,runame) for i in 1:T, k in 1:nsteps+1, f in 1:Nf]
    pptilde_t_mc = real.(sum(pptilde_t_mc,dims = 2)[:,1,:,:]./Nsff)
    global pptilde_t = [uwreal(pptilde_t_mc[i,k,f] ,runame) for i in 1:T, k in 1:nsteps+1, f in 1:Nf]

    pp_corr_tfl_mc = sum(pp_corr_tfl_mc,dims = 2)[:,1,:,:]./Nsff
    global pp_corr_tfl = [uwreal(pp_corr_tfl_mc[i,k,f] ,runame) for i in 1:T, k in 1:Nfl, f in 1:Nf]
    ap_corr_tfl_mc = real.(sum(ap_corr_tfl_mc,dims = 2)[:,1,:,:]./Nsff)
    global ap_corr_tfl = [uwreal(ap_corr_tfl_mc[i,k,f] ,runame) for i in 1:T, k in 1:Nfl, f in 1:Nf]

    Sigma_mc = real.(sum(Sigma_mc,dims = 2)[:,1,:,:]./Nsbf)
    global Sigma = [uwreal(Sigma_mc[i,k,f] ,runame) for i in 1:T, k in 1:Nfl, f in 1:Nf]
    Sigma_cfl_mc = real.(sum(Sigma_cfl_mc,dims = 2)[:,1,:,:]./Nsbf)
    global Sigma_cfl = [uwreal(Sigma_cfl_mc[i,k,f] ,runame) for i in 1:T, k in 1:Nfl, f in 1:Nf]

    Sigma2_mc = real.(sum(Sigma2_mc,dims = 2)[:,1,:,:]./Nsbf)
    global Sigma2 = [uwreal(Sigma2_mc[i,k,f] ,runame) for i in 1:T, k in 1:Nfl, f in 1:Nf]
    Sigma2_cfl_mc = real.(sum(Sigma2_cfl_mc,dims = 2)[:,1,:,:]./Nsbf)
    global Sigma2_cfl = [uwreal(Sigma2_cfl_mc[i,k,f] ,runame) for i in 1:T, k in 1:Nfl, f in 1:Nf]

    ChiDchi_mc = real.(sum(ChiDchi_mc,dims = 2)[:,1,:,:]./Nsbf)
    global ChiDchi = [uwreal(ChiDchi_mc[i,k,f] ,runame) for i in 1:T, k in 1:Nfl, f in 1:Nf]
    ChiDchi_cfl_mc = real.(sum(ChiDchi_cfl_mc,dims = 2)[:,1,:,:]./Nsbf)
    global ChiDchi_cfl = [uwreal(ChiDchi_cfl_mc[i,k,f] ,runame) for i in 1:T, k in 1:Nfl, f in 1:Nf]

    return nothing
end