sfcf.jl/src/io.jl

using LatticeGPU
using TOML
using ArgParse
using CUDA
using BDIO
using InteractiveUtils
using MD5


"""
    function read_input()

Stores as global variables 'parsed_args' (info from the command line) and 'params' (info from the input file)
"""
function read_input()
    global parsed_args = parse_commandline()
    global params = TOML.parsefile(parsed_args["i"])
    return nothing
end


function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "-i"
        help = "Input parameters file"
        required = true
        arg_type = String

        "-c"
        help = "Gauge configuration file"
        required = true
        arg_type = String

        "--cern"
        help = "Config written with the export_cnfg_cern() convention"
        action = :store_true

    end

    return parse_args(s)
end


"""
    function load_gauge_field()

Returns the gauge field and computes the Csw term
"""
function load_gauge_field()
    if !parsed_args["cern"]
    println(log_file,"\nReading gauge field with LGPU native convention from: ", parsed_args["c"], "...")
        U,_ = read_cnfg(parsed_args["c"])
    else
    println(log_file,"\nReading gauge field with CERN convention from: ", parsed_args["c"], "...")
        U = read_cnfg_cern(parsed_args["c"],lp)
    end

    f = open(parsed_args["c"],"r")
    println(log_file,"MD5 checksum of gauge config: ",bytes2hex(md5(f)),"\n")
    close(f)

    Csw!(dws, U, gp, lp)
    return U
end

function read_cnfg_cern(path::String,lp::SpaceParm)
    Ucpu = Array{SU3{Float64}, 3}(undef, lp.bsz, lp.ndim, lp.rsz)
    file = open(path)
    for t in 1:lp.iL[4]
        for i in 1:lp.iL[1]
            for j in 1:lp.iL[2]
                for k in 1:lp.iL[3]

                    for d in [4,1,2,3]

                        f,r = point_index(CartesianIndex((i,j,k,t)),lp)

                        #a11 !!

                        re11 = read(file,Float64)
                        co11 = read(file,Float64)

                        #a12 !!

                        re12 = read(file,Float64)
                        co12 = read(file,Float64)

                        #a13 !!

                        re13 = read(file,Float64)
                        co13 = read(file,Float64)

                        #a21 !!

                        re21 = read(file,Float64)
                        co21 = read(file,Float64)

                        #a22 !!

                        re22 = read(file,Float64)
                        co22 = read(file,Float64)

                        #a23 !!

                        re23 = read(file,Float64)
                        co23 = read(file,Float64)

                        #a31

                        re31 = read(file,Float64)
                        co31 = read(file,Float64)

                        #a32

                        re32 = read(file,Float64)
                        co32 = read(file,Float64)

                        #a33

                        re33 = read(file,Float64)
                        co33 = read(file,Float64)

                        (Ucpu[f,d,r] = SU3{Float64}(re11 + im*co11, re12 + im*co12, re13 + im*co13,
                                                    re21 + im*co21, re22 + im*co22, re23 + im*co23))
                    end
                end
            end
        end
    end
    length(read(file)) == (prod(lp.iL[1:3])*4*8*9*2) ? nothing : error("File not fully read")
    close(file)

    return CuArray(Ucpu)
end

"""
    function load_structs()

Stores in global variables the needed structures, i.e. lp, gp, dpar, dws, ymws
"""
function load_structs()

    global lp = SpaceParm{4}(tuple(params["Space"]["size"]...), tuple(params["Space"]["blocks"]...),BC_SF_ORBI, (0,0,0,0,0,0))
    global gp = GaugeParm{Float64}(SU3{Float64},params["Fermion"]["beta"],1.0,(params["Space"]["cG"],0.0),params["Space"]["phiT"],lp.iL);
    global dpar = DiracParam{Float64}(SU3fund,(1/(2*params["Fermion"]["kappa"])) - 4,params["Fermion"]["csw"],ntuple(i -> exp((i!=4)*im*params["Fermion"]["theta"]/lp.iL[i]),4),0.0,params["Fermion"]["ct"]);
    global dws = DiracWorkspace(SU3fund,Float64,lp);
    global ymws = YMworkspace(SU3,Float64,lp);

    return nothing
end

function write_log()

    println(log_file,"Running sfcf.jl by ", params["Run"]["user"],". Name of the run: ",params["Run"]["name"])
    println(log_file,"")

    print(log_file,"Calling: ")
    print(log_file,PROGRAM_FILE*" "); for x in ARGS; print(log_file,x*" "); end
    println(log_file,"\n")

    println(log_file,"Version info:")
    versioninfo(log_file)
    println(log_file,"")

    println(log_file,"Reading input file from:", parsed_args["i"], "...\n")

    println(log_file,"Parameters:")
    println(log_file,"Lattice size:  ", lp.iL)
    println(log_file,"Phi0         = ", params["Space"]["phi0"])
    println(log_file,"PhiT         = ", params["Space"]["phiT"])
    println(log_file,"cG           = ", gp.cG[1])
    println(log_file,"kappa        = ", params["Fermion"]["kappa"])
    println(log_file,"theta        = ", params["Fermion"]["theta"])
    println(log_file,"csw          = ", dpar.csw)
    println(log_file,"ct           = ", dpar.ct)
    println(log_file,"tolerance    = ", params["Solver"]["tolerance"])
    println(log_file,"maxiter      = ", params["Solver"]["maxiter"])

    flush(log_file)
    return nothing
end


function save_correlators()

    ihdr = [convert(Int32,1708683512)]
    fname = "./output/"*params["Run"]["name"]*".bdio"

    if isfile(fname)
        fb = BDIO_open(fname, "a")
        println(log_file,"\n\nAppending output to "*fname*"\n\n")
    else
        fb = BDIO_open(fname, "w", "BDIO output from sfcf.jl")
        println(log_file,"Creating new BDIO output file "*fname)

        BDIO_start_record!(fb, BDIO_BIN_GENERIC, 14)
        BDIO_write!(fb, ihdr)
        BDIO_write_hash!(fb)

        BDIO_start_record!(fb, BDIO_BIN_GENERIC, 1)
        BDIO_write!(fb, [convert(Int32, 4)])
        BDIO_write!(fb, [convert(Int32, lp.iL[i]) for i in 1:4])
        BDIO_write!(fb, [convert(Int32, lp.ntw[i]) for i in 1:6])
        BDIO_write!(fb, [dpar.m0, dpar.csw, dpar.tm, dpar.ct])
        BDIO_write!(fb, [dpar.th[i] for i in 1:4])

        BDIO_write_hash!(fb)
    end


    BDIO_start_record!(fb, BDIO_BIN_GENERIC, 8)
    BDIO_write!(fb,basename(parsed_args["c"]))
    fg = open(parsed_args["c"],"r")
    BDIO_write!(fb,bytes2hex(md5(fg)))
    close(fg)

    BDIO_write!(fb,fP)
    BDIO_write!(fb,fA)
    BDIO_write!(fb,[f1])
    BDIO_write!(fb,gP)
    BDIO_write!(fb,gA)
    BDIO_write!(fb,kV)
    BDIO_write!(fb,lV)
    BDIO_write!(fb,[k1])
    BDIO_write!(fb,kT)
    BDIO_write!(fb,lT)

    BDIO_write_hash!(fb)

    BDIO_close!(fb)

    return nothing
end