Write I/O functions based on ccall to cut allocations

MetaAnalysis/src/MetaAnalysis.jl

@@ -18,7 +18,13 @@ export Bootstrap, Jackknife, UWerr, analyze, modify_bias
 
 abstract type AbstractErrorEstimator end
 
+<<<<<<< Updated upstream
 include("bias.jl")
+=======
+# TODO: get rid of MetaQCD dependency
+# include("../../src/bias/metadynamics.jl")
+# include("../../src/bias/opes.jl")
+>>>>>>> Stashed changes
 include("viz.jl")
 include("autocorr.jl")
 include("bootstrap.jl")

MetaAnalysis/src/biasmod.jl

@@ -1,8 +1,13 @@
+<<<<<<< Updated upstream
 function modify_bias(b::MetaBias, order; savefile="", cvlims=b.bias.cvlims, L=80)
+=======
+function modify_bias(b::MetaBias; order=1, cvlims=b.bias.cvlims, L=80)
+>>>>>>> Stashed changes
     x = cvlims[1]:0.01:cvlims[2]
     yn = b.bias.(x)
     yt, ys = SingularSpectrumAnalysis.analyze(yn, L, robust=true)
     # TODO: Fit and return parameters
+<<<<<<< Updated upstream
     p0 = zeros(2order+1)
     p0 .= 1.0
     p0[1] = 5.0
@@ -25,10 +30,19 @@ function modify_bias(b::MetaBias, order; savefile="", cvlims=b.bias.cvlims, L=80
 
     model = barrier_func(order)
     bfit = curve_fit(model, x, y, p0)
+=======
+    p0 = zeros(order)
+    p0 .= 1
+    y = ys[:, 1] .- minimum(view(ys, :, 1))
+    model = fitting_func(order)
+    bfit = curve_fit(model, x, y, p0)
+    # display(plot(x, yt))
+>>>>>>> Stashed changes
     plt = plot(x, y)
     fity = model(x, coef(bfit))
     plot!(plt, x, fity)
     display(plt)
+<<<<<<< Updated upstream
     return (yt, ys), coef(bfit), sum(bfit.resid.^2)
 end
 
@@ -40,5 +54,17 @@ function barrier_func(order)
             out .+= p[2(i-1)+2] * cospi.(2p[2(i-1)+3]*x)
         end
         return out
+=======
+    return coef(bfit), sum(bfit.resid.^2)
+end
+
+function fitting_func(order)
+    @assert order ≥ 1
+    function f(x, p)
+        out = p[1] * cospi(p[2]*x).^2
+        for i in 2:order
+            out .+= p[2(i-1)+1] * cospi.(p[2(i-1)+2]*x).^2
+        end
+>>>>>>> Stashed changes
     end
 end

Project.toml

@@ -10,13 +10,13 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ArnoldiMethod = "ec485272-7323-5ecc-a04f-4719b315124d"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+Format = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
-MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
@@ -26,6 +26,7 @@ RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+StaticTools = "86c06d3c-3f03-46de-9781-57580aa96d0a"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StrideArraysCore = "7792a7ef-975c-4747-a70f-980b88e8d1da"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

metaqcd_sim.jl

@@ -5,7 +5,8 @@ using MetaQCD.MPI
 using MetaQCD: run_sim
 using MetaQCD.Parameters: lower_case
 
-@assert MPI.Comm_size(MPI.COMM_WORLD) == 1 "metaqcd_sim can not be used with mpi, only metaqcd_build"
+const COMM = MPI.COMM_WORLD
+@assert MPI.Comm_size(COMM) == 1 "metaqcd_sim can not be used with mpi, only metaqcd_build"
 
 if length(ARGS) == 0
     error("""

src/MetaQCD.jl

@@ -52,7 +52,7 @@ import .Fields: normalize!, plaquette, plaquette_trace_sum, random_gauges!
 import .Fields: staple, staple_eachsite!, wilsonloop, to_backend
 import .Fields: Tensorfield, calc_kinetic_energy, gaussian_TA!
 import .Fields: Fermionfield, gaussian_pseudofermions!
-import .Measurements: measure, get_value, top_charge
+import .Measurements: measure, top_charge
 import .Measurements: EnergyDensityMeasurement, GaugeActionMeasurement, PlaquetteMeasurement
 import .Measurements: PolyakovMeasurement, TopologicalChargeMeasurement
 import .Measurements: WilsonLoopMeasurement
@@ -78,7 +78,7 @@ export StaggeredEOPreDiracOperator, even_odd, sample_pseudofermions!
 export StaggeredFermionAction, StaggeredEOPreFermionAction
 export WilsonFermionAction, WilsonEOPreFermionAction
 export calc_fermion_action, gaussian_pseudofermions!
-export measure, get_value, top_charge
+export measure, top_charge
 export EnergyDensityMeasurement, GaugeActionMeasurement, PlaquetteMeasurement
 export PolyakovMeasurement, TopologicalChargeMeasurement, WilsonLoopMeasurement
 export ParameterSet, construct_params_from_toml

src/bias/Bias.jl

@@ -37,10 +37,9 @@ struct Bias{TCV,TS,TB,TW,T}
     is_static::Bool
     bias::TB
     kinds_of_weights::TW
-    biasfile::String
-    datafile::String
+    biasfile::T
+    datafile::T
     write_bias_every::Int64
-    fp::T
 end
 
 function Bias(p::ParameterSet, U; use_mpi=false, instance=1)
@@ -69,34 +68,28 @@ function Bias(p::ParameterSet, U; use_mpi=false, instance=1)
         is_opes = bias isa OPES
         kinds_of_weights = is_opes ? ["opes"] : p.kinds_of_weights
         ext = is_opes ? "opes" : "metad"
-        biasfile = MYRANK==0 ? p.bias_dir * "/stream_$(inum).$(ext)" : ""
-        datafile = p.measure_dir * "/bias_data_$inum.txt"
+        biasfile = MYRANK==0 ? joinpath(p.bias_dir, "stream_$(inum).$(ext)") : ""
+        datafile = joinpath(p.measure_dir, "bias_data_$inum.txt")
         # FIXME: For some reason this errors with MPI on the UNI's cluster
-        fp = open(datafile, "w")
-        str = @sprintf("%-9s\t%-22s", "itrj", "cv")
-        print(fp, str)
+        open(datafile, "w") do fp
+            @printf(fp, "%-11s%-25s", "itrj", "cv")
 
-        for name in kinds_of_weights
-            str = @sprintf("\t%-22s", "weight_$(name)")
-            print(fp, str)
+            for name in kinds_of_weights
+                @printf(fp, "%-25s", "weight_$(name)")
+            end
+            cnewline(fp)
         end
-
-        println(fp)
     elseif bias isa Parametric
         kinds_of_weights = ["branduardi"]
         biasfile = ""
-        datafile = p.measure_dir * "/bias_data_$inum.txt"
-        fp = open(datafile, "w")
-        str = @sprintf("%-9s\t%-22s\t%-22s\n", "itrj", "cv", "weight_branduardi")
-        println(fp, str)
+        datafile = joinpath(p.measure_dir, "bias_data_$inum.txt")
+        open(datafile, "w") do fp
+            @printf(fp, "%-11s%-25s%-25s", "itrj", "cv", "weight_branduardi")
+            cnewline(fp)
+        end
         @level1(
             "|  @info: Parametric bias defaults to static and weight-type \"branduardi\""
         )
-    else
-        biasfile = ""
-        datafile = ""
-        kinds_of_weights = nothing
-        fp = nothing
     end
 
     @level1("|  BIASFILE: $(biasfile)")
@@ -109,7 +102,8 @@ function Bias(p::ParameterSet, U; use_mpi=false, instance=1)
 
     # write to file after construction to make sure nothing went wrong
     MYRANK == 0 && write_to_file(bias, biasfile)
-    @level1("└\n")
+    @level1("└")
+    @level1("")
     return Bias(
         TCV(),
         smearing,
@@ -119,7 +113,6 @@ function Bias(p::ParameterSet, U; use_mpi=false, instance=1)
         biasfile,
         datafile,
         write_bias_every,
-        fp,
     )
 end
 
@@ -140,8 +133,6 @@ end
 
 kind_of_cv(::NoBias) = nothing
 kind_of_cv(b::Bias) = b.kind_of_cv
-Base.close(::NoBias) = nothing
-Base.close(b::Bias{TCV,TS,TB,T}) where {TCV,TS,TB,T} = T ≢ Nothing ? close(b.fp) : nothing
 update_bias!(::NoBias, args...) = nothing
 update_bias!(::Nothing, args...) = nothing
 write_to_file(::AbstractBias, args...) = nothing

src/bias/opes.jl

@@ -523,7 +523,7 @@ const kernel_header = "#$(rpad("height", 20))\t$(rpad("center", 20))\t$(rpad("si
 write_to_file(::OPES, ::Nothing) = nothing
 
 function write_to_file(o::OPES, filename::String)
-    filename == "" && return nothing
+    isnothing(filename) && return nothing
     (tmppath, tmpio) = mktemp()
     print(tmpio, "#")
     [print(tmpio, "$(var)\t") for var in state_vars]

src/bias/parametric.jl

@@ -42,12 +42,12 @@ function (p::Parametric)(cv)
     Q′ = Z * cv
 
     if in_bounds(cv, lb, ub)
-        return Q * cv^2 + A * cos(2π * Q′)
+        return Q * cv^2 + A * cospi(Q′)^2
     elseif cv < lb
-        penalty = Q * cv^2 + A * cos(2π * Q′) + pen * (cv - lb)^2
+        penalty = Q * cv^2 + A * cospi(Q′)^2 + pen * (cv - lb)^2
         return penalty
     else
-        penalty = Q * cv^2 + A * cos(2π * Q′) + pen * (cv - ub)^2
+        penalty = Q * cv^2 + A * cospi(Q′)^2 + pen * (cv - ub)^2
         return penalty
     end
 end
@@ -59,19 +59,19 @@ function ∂V∂Q(p::Parametric, cv)
     Q′ = Z * cv
 
     if in_bounds(cv, lb, ub)
-        return 2Q * cv - 2π * Z * A * sin(2π * Q′)
+        return 2Q * cv - 2π * Z * A * sinpi(Q′)
     elseif cv < lb
-        penalty = 2Q * cv - 2π * Z * A * sin(2π * Q′) + 2pen * (cv - lb)
+        penalty = 2Q * cv - 2π * Z * A * sinpi(Q′) + 2pen * (cv - lb)
         return penalty
     else
-        penalty = 2Q * cv - 2π * Z * A * sin(2π * Q′) + 2pen * (cv - ub)
+        penalty = 2Q * cv - 2π * Z * A * sinpi(Q′) + 2pen * (cv - ub)
         return penalty
     end
 end
 
 function integral(p::Parametric, lb, ub) # XXX: Why does this function exist?
     Q, A, Z = p.Q, p.A, p.Z
-    num = 3A * sin(2π*Z*ub) + 2π*Q*Z*ub^3 - 3A * sin(2π*Z*lb) - 2π*Q*Z*lb^3
+    num = 3A * sinpi(Z*ub) + 2π*Q*Z*ub^3 - 3A * sinpi(Z*lb) - 2π*Q*Z*lb^3
     denom = 6π * Z
     return num / denom
 end

src/bias/weights.jl

@@ -5,27 +5,28 @@ https://pubs.acs.org/doi/pdf/10.1021/acs.jctc.9b00867
 calc_weights(::Nothing, args...) = nothing
 calc_weights(::NoBias, args...) = nothing
 
-function calc_weights(fps, b::Vector{<:Bias}, cv, itrj)
+function calc_weights(filenames, b::Vector{<:Bias}, cv, itrj)
     for i in eachindex(b)
-        calc_weights(fps[i], b[i], cv[i], itrj)
+        calc_weights(filenames[i], b[i], cv[i], itrj)
     end
 
     return nothing
 end
 
-function calc_weights(fp, b::Bias{TCV,TS,TB}, cv, itrj) where {TCV,TS,TB}
-    b.kinds_of_weights === nothing && return nothing
-    str = @sprintf("%-9i\t%+-22.15E", itrj, cv)
-    for method in b.kinds_of_weights
-        w = calc_weight(b.bias, cv, method)
-        str *= @sprintf("\t%-22.15E", w)
-    end
-
-    @level1(str * " # cv weight")
+function calc_weights(filename, b::Bias{TCV,TS,TB}, cv, itrj) where {TCV,TS,TB}
+    if isnothing(filename)
+        for method in b.kinds_of_weights
+            w = calc_weight(b.bias, cv, method)
+            @level1("$itrj\t$cv\t$w # cv weight_$method")
+        end
+    else
+        fp = copen(filename, "a")
+        cprint(fp, "%-11i%+-25.15E", itrj, cv)
 
-    if fp ≢ nothing
-        println(fp, str)
-        flush(fp)
+        for method in b.kinds_of_weights
+            w = calc_weight(b.bias, cv, method)
+            cprint(fp, "%-25.15E", w)
+        end
     end
 
     return nothing

src/diracoperators/staggered.jl

@@ -186,8 +186,8 @@ function calc_fermion_action(
     n = get_n(rhmc)
     D = fermion_action.D(U)
     DdagD = DdaggerD(D)
-    ψs = fermion_action.rhmc_temps1[1:n+1]
-    ps = fermion_action.rhmc_temps2[1:n+1]
+    ψs = fermion_action.rhmc_temps1
+    ps = fermion_action.rhmc_temps2
     temp1, temp2 = fermion_action.cg_temps
 
     for v in ψs
@@ -225,8 +225,8 @@ function sample_pseudofermions!(ϕ, fermion_action::StaggeredFermionAction{Nf},
     n = get_n(rhmc)
     D = fermion_action.D(U)
     DdagD = DdaggerD(D)
-    ψs = fermion_action.rhmc_temps1[1:n+1]
-    ps = fermion_action.rhmc_temps2[1:n+1]
+    ψs = fermion_action.rhmc_temps1
+    ps = fermion_action.rhmc_temps2
     temp1, temp2 = fermion_action.cg_temps
 
     for v in ψs

src/diracoperators/wilson.jl

@@ -212,8 +212,8 @@ function calc_fermion_action(
     n = get_n(rhmc)
     D = fermion_action.D(U)
     DdagD = DdaggerD(D)
-    ψs = fermion_action.rhmc_temps1[1:n+1]
-    ps = fermion_action.rhmc_temps2[1:n+1]
+    ψs = fermion_action.rhmc_temps1
+    ps = fermion_action.rhmc_temps2
     temp1, temp2 = fermion_action.cg_temps
 
     for v in ψs
@@ -251,8 +251,8 @@ function sample_pseudofermions!(ϕ, fermion_action::WilsonFermionAction{Nf}, U)
     n = get_n(rhmc)
     D = fermion_action.D(U)
     DdagD = DdaggerD(D)
-    ψs = fermion_action.rhmc_temps1[1:n+1]
-    ps = fermion_action.rhmc_temps2[1:n+1]
+    ψs = fermion_action.rhmc_temps1
+    ps = fermion_action.rhmc_temps2
     temp1, temp2 = fermion_action.cg_temps
 
     for v in ψs

src/diracoperators/wilson_eo.jl

@@ -258,7 +258,7 @@ function calc_fermion_action(
 
     clear!(ψ_eo) # initial guess is zero
     solve_dirac!(ψ_eo, DdagD, ϕ_eo, temp1, temp2, temp3, cg_tol, cg_maxiters) # ψ = (D†D)⁻¹ϕ
-    Sf = dot(ϕ_eo, ψ_eo) #= - 2trlog(D.D_diag) =#
+    Sf = #= dot(ϕ_eo, ψ_eo) =# - 2trlog(D.D_diag)
     return real(Sf)
 end
 

src/fields/gpu_kernels/utils.jl

@@ -59,8 +59,6 @@ end
     return val
 end
 
-using MacroTools: @capture, prewalk, postwalk, splitdef, combinedef
-
 # @inline function tile(NX::Integer, NY::Integer, NZ::Integer, NT::Integer)
 #     maxn = 1024
 #     num1 = min(maxn, NX)

src/forces/forces.jl

@@ -1,3 +1,7 @@
+# XXX: Maybe make this its own module?
+
+# module Forces
+
 import ..DiracOperators: StaggeredDiracOperator, StaggeredFermionAction
 import ..DiracOperators: StaggeredEOPreDiracOperator, StaggeredEOPreFermionAction
 import ..DiracOperators: WilsonDiracOperator, WilsonFermionAction
@@ -45,3 +49,5 @@ include("gpu_kernels/gauge_force.jl")
 include("gpu_kernels/bias_force.jl")
 include("gpu_kernels/wilson_force.jl")
 include("gpu_kernels/staggered_force.jl")
+
+# end

src/forces/staggered_eo_force.jl

@@ -30,8 +30,8 @@ function calc_dSfdU!(
     D = fermion_action.D(U)
     DdagD = DdaggerD(D)
     anti = D.anti_periodic
-    Xs = fermion_action.rhmc_temps1[1:n+1]
-    Ys = fermion_action.rhmc_temps2[1:n+1]
+    Xs = fermion_action.rhmc_temps1
+    Ys = fermion_action.rhmc_temps2
     temp1, temp2 = fermion_action.cg_temps
 
     for X in Xs

src/forces/staggered_force.jl

@@ -29,8 +29,8 @@ function calc_dSfdU!(
     D = fermion_action.D(U)
     DdagD = DdaggerD(D)
     anti = D.anti_periodic
-    Xs = fermion_action.rhmc_temps1[1:n+1]
-    Ys = fermion_action.rhmc_temps2[1:n+1]
+    Xs = fermion_action.rhmc_temps1
+    Ys = fermion_action.rhmc_temps2
     temp1, temp2 = fermion_action.cg_temps
 
     for X in Xs