<commit_before>## staghmc_s.nim (hypsmear) with Hasenbusch masses.

import qex, gauge, gauge/hypsmear, physics/qcdTypes, physics/stagSolve
import mdevolve
import math, sequtils, strutils, times

const ReversibilityCheck {.booldefine.} = false

type IntProc = proc(T,V:Integrator; steps:int):Integrator
converter toIntProc(s:string):IntProc =
  template mkProc1(s:untyped):IntProc =
    proc mkInt(T,V:Integrator; steps:int):Integrator {.gensym.} =
      `mk s`(T = T, V = V, steps = steps)
    mkInt
  template mkProc2(s:untyped):IntProc =
    proc mkInt(T,V:Integrator; steps:int):Integrator {.gensym.} =
      `mk s`(T = T, V = V, steps = steps, ss[1].parseFloat)
    mkInt
  template mkProc3(s:untyped):IntProc =
    proc mkInt(T,V:Integrator; steps:int):Integrator {.gensym.} =
      `mk s`(T = T, V = V, steps = steps, ss[1].parseFloat, ss[2].parseFloat)
    mkInt
  template mkProc4(s:untyped):IntProc =
    proc mkInt(T,V:Integrator; steps:int):Integrator {.gensym.} =
      `mk s`(T = T, V = V, steps = steps, ss[1].parseFloat, ss[2].parseFloat, ss[3].parseFloat)
    mkInt
  template mkProc5(s:untyped):IntProc =
    proc mkInt(T,V:Integrator; steps:int):Integrator {.gensym.} =
      `mk s`(T = T, V = V, steps = steps, ss[1].parseFloat, ss[2].parseFloat, ss[3].parseFloat, ss[4].parseFloat)
    mkInt
  let ss = s.split(',')
  # Omelyan's triple star integrators, see Omelyan et. al. (2003)
  case ss[0]:
  of "2MN":
    if ss.len == 1: return mkProc1(Omelyan2MN)
    else: return mkProc2(Omelyan2MN)
  of "4MN5FP":
    if ss.len == 1: return mkProc1(Omelyan4MN5FP)
    elif ss.len == 2: return mkProc2(Omelyan4MN5FP)
    elif ss.len == 3: return mkProc3(Omelyan4MN5FP)
    elif ss.len == 4: return mkProc4(Omelyan4MN5FP)
    elif ss.len == 5: return mkProc5(Omelyan4MN5FP)
    else: return mkProc2(Omelyan4MN5FP)
  of "4MN5FV":
    if ss.len == 1: return mkProc1(Omelyan4MN5FV)
    elif ss.len == 2: return mkProc2(Omelyan4MN5FV)
    elif ss.len == 3: return mkProc3(Omelyan4MN5FV)
    elif ss.len == 4: return mkProc4(Omelyan4MN5FV)
    elif ss.len == 5: return mkProc5(Omelyan4MN5FV)
    else: return mkProc2(Omelyan4MN5FV)
  of "6MN7FV": return mkProc1(Omelyan6MN7FV)
  of "4MN3F1GP":  # lambda = 0.2725431326761773  is  FUEL f3g a0=0.109
    if ss.len == 1: return mkProc1(Omelyan4MN3F1GP)
    else: return mkProc2(Omelyan4MN3F1GP)
  of "4MN4F2GVG": return mkProc1(Omelyan4MN4F2GVG)
  of "4MN4F2GV": return mkProc1(Omelyan4MN4F2GV)
  of "4MN5F1GV": return mkProc1(Omelyan4MN5F1GV)
  of "4MN5F1GP": return mkProc1(Omelyan4MN5F1GP)
  of "4MN5F2GV": return mkProc1(Omelyan4MN5F2GV)
  of "4MN5F2GP": return mkProc1(Omelyan4MN5F2GP)
  of "6MN5F3GP": return mkProc1(Omelyan6MN5F3GP)
  else:
    qexError "Cannot parse integrator: '", s, "'\n",
      """Available integrators (with default parameters):
      2MN,0.1931833275037836
      4MN5FP,0.2750081212332419,-0.1347950099106792,-0.08442961950707149,0.3549000571574260
      4MN5FV,0.2539785108410595,-0.03230286765269967,0.08398315262876693,0.6822365335719091
      6MN7FV
      4MN3F1GP,0.2470939580390842
      4MN4F2GVG
      4MN4F2GV
      4MN5F1GV
      4MN5F1GP
      4MN5F2GV
      4MN5F2GP
      6MN5F3GP"""

qexinit()

letParam:
  lat = @[8,8,8,8]
  beta = 6.0
  adjFac = -0.25
  tau = 2.0
  trajs = 10
  seed:uint64 = int(1000*epochTime())
  gintalg:IntProc = "4MN5F2GP"
  gsteps = 4
  mass = @[0.1]  # mass for each staggered species
  hmasses0 = @[0.2,0.4]  # Hasenbusch masses for mass[0]
  hmasses1 = if mass.len>1: hmasses0 else: @[]  # Hasenbusch masses for mass[1]
  hmasses2 = if mass.len>2: hmasses0 else: @[]  # Hasenbusch masses for mass[2]
  hmasses3 = if mass.len>3: hmasses0 else: @[]  # Hasenbusch masses for mass[3]
  hmasses4 = if mass.len>4: hmasses0 else: @[]  # Hasenbusch masses for mass[4]
  fintalg:IntProc = "4MN5F2GP"
  fsteps = repeat(4, mass.len)  # nsteps for each mass
  hfsteps0 = fsteps[0].repeat hmasses0.len  # nsteps for Hasenbusch masses 0
  hfsteps1 = fsteps[0].repeat hmasses1.len  # nsteps for Hasenbusch masses 1
  hfsteps2 = fsteps[0].repeat hmasses2.len  # nsteps for Hasenbusch masses 2
  hfsteps3 = fsteps[0].repeat hmasses3.len  # nsteps for Hasenbusch masses 3
  hfsteps4 = fsteps[0].repeat hmasses4.len  # nsteps for Hasenbusch masses 4
  arsq = 1e-20  # CG r^2 for fermion action
  frsq = repeat(1e-12, mass.len)  # CG r^2 for fermion force for each mass
  hfrsq0 = frsq[0].repeat hmasses0.len  # frsq for Hasenbusch masses 0
  hfrsq1 = frsq[0].repeat hmasses1.len  # frsq for Hasenbusch masses 1
  hfrsq2 = frsq[0].repeat hmasses2.len  # frsq for Hasenbusch masses 2
  hfrsq3 = frsq[0].repeat hmasses3.len  # frsq for Hasenbusch masses 3
  hfrsq4 = frsq[0].repeat hmasses4.len  # frsq for Hasenbusch masses 4
  pbpmass = mass
  pbpreps = repeat(1, pbpmass.len)
  pbprsq = arsq
  maxits = 10000
  useFG2:bool = 0

echoParams()

if mass.len > 5:
  qexError "Unimlemented for mass: ", mass
if mass.len != fsteps.len or
    mass.len != frsq.len:
  qexError "Parameters for staggered species mismatch."

let
  hmasses = @[hmasses0, hmasses1, hmasses2, hmasses3, hmasses4][0..<mass.len]
  hfsteps = @[hfsteps0, hfsteps1, hfsteps2, hfsteps3, hfsteps4][0..<mass.len]
  hfrsq = @[hfrsq0, hfrsq1, hfrsq2, hfrsq3, hfrsq4][0..<mass.len]

for k in 0..<mass.len:
  if hmasses[k].len != hfsteps[k].len or
      hmasses[k].len != hfrsq[k].len:
    qexError "Hasenbusch parameters lengths mismatch."

if pbpmass.len != pbpreps.len:
  qexError "The lengths of pbpmass and pbpreps differ."

let
  lo = lat.newLayout
  #gc = GaugeActionCoeffs(plaq:6)
  gc = GaugeActionCoeffs(plaq: beta, adjplaq: beta*adjFac)
  vol = lo.physVol

var r = lo.newRNGField(RngMilc6, seed)
var R:RngMilc6  # global RNG
R.seed(seed, 987654321)

var g = lo.newgauge
#g.random r
g.unit

echo 6.0*g.plaq
echo gc.actionA g

var
  p = lo.newgauge
  f = lo.newgauge
  g0 = lo.newgauge
  sg0 = lo.newgauge

var ftmp = lo.ColorVector()
type CV = typeof(ftmp)
var
  phi = newseq[seq[CV]](mass.len)
  psi = newseq[seq[CV]](mass.len)
for k in 0..<mass.len:
  phi[k] = newseq[typeof(ftmp)](hmasses[k].len+1)
  psi[k] = newseq[typeof(ftmp)](hmasses[k].len+1)
  for i in 0..<phi[k].len:
    phi[k][i] = lo.ColorVector()
    psi[k][i] = lo.ColorVector()

var pbpsp = initSolverParams()
pbpsp.r2req = pbprsq
pbpsp.maxits = maxits

var spa = initSolverParams()
#spa.subsetName = "even"
spa.r2req = arsq
spa.maxits = maxits

var spf = newseq[typeof(spa)](mass.len)
for k in 0..<mass.len:
  spf[k] = initSolverParams()
  #spf[k].subsetName = "even"
  spf[k].r2req = frsq[k]
  spf[k].maxits = maxits
  spf[k].verbosity = 0

var spfh = newseq[typeof(spf)](mass.len)
for k in 0..<mass.len:
  spfh[k] = newseq[typeof(spa)](hfrsq[k].len)
  for i in 0..<spfh[k].len:
    spfh[k][i] = initSolverParams()
    spfh[k][i].r2req = hfrsq[k][i]
    spfh[k][i].maxits = maxits
    spfh[k][i].verbosity = 0

var
  info: PerfInfo
  coef = HypCoefs(alpha1:0.4, alpha2:0.5, alpha3:0.5)
echo "smear = ",coef

var
  sg = lo.newGauge  # smeared gauge
  gg = lo.newgauge  # FG backup gauge
  sgg = lo.newgauge  # FG backup smeared gauge
let
  stag = newStag(sg)
  stag0 = newStag(sg0)
  stagg = newStag(sgg)

proc pbp(stag:any) =
  var ftmp2 = lo.ColorVector()
  for k in 0..<pbpmass.len:
    for i in 0..<pbpreps[k]:
      threads:
        ftmp.u1 r
      stag0.solve(ftmp2, ftmp, pbpmass[k], pbpsp)
      threads:
        var pbp = ftmp2.norm2
        threadMaster:
          echo "MEASpbp mass ",pbpmass[k]," : ",0.5*pbpmass[k]*pbp/vol.float

proc fgsave =
  threads:
    for mu in 0..<g.len:
      gg[mu] := g[mu]
proc fgload =
  threads:
    for mu in 0..<g.len:
      g[mu] := gg[mu]

proc smearRephase(g: any, sg: any):auto {.discardable.} =
  tic()
  let smearedForce = coef.smear(g, sg, info)
  toc("smear")
  threads:
    sg.setBC
    threadBarrier()
    sg.stagPhase
  toc("BC & Phase")
  smearedForce

template pnorm2(p2:float) =
  threads:
    var p2t = 0.0<commit_msg>staghmc_sh: fix pbp measurement<commit_after>## staghmc_s.nim (hypsmear) with Hasenbusch masses.

import qex, gauge, gauge/hypsmear, physics/qcdTypes, physics/stagSolve
import mdevolve
import math, sequtils, strutils, times

const ReversibilityCheck {.booldefine.} = false

type IntProc = proc(T,V:Integrator; steps:int):Integrator
converter toIntProc(s:string):IntProc =
  template mkProc1(s:untyped):IntProc =
    proc mkInt(T,V:Integrator; steps:int):Integrator {.gensym.} =
      `mk s`(T = T, V = V, steps = steps)
    mkInt
  template mkProc2(s:untyped):IntProc =
    proc mkInt(T,V:Integrator; steps:int):Integrator {.gensym.} =
      `mk s`(T = T, V = V, steps = steps, ss[1].parseFloat)
    mkInt
  template mkProc3(s:untyped):IntProc =
    proc mkInt(T,V:Integrator; steps:int):Integrator {.gensym.} =
      `mk s`(T = T, V = V, steps = steps, ss[1].parseFloat, ss[2].parseFloat)
    mkInt
  template mkProc4(s:untyped):IntProc =
    proc mkInt(T,V:Integrator; steps:int):Integrator {.gensym.} =
      `mk s`(T = T, V = V, steps = steps, ss[1].parseFloat, ss[2].parseFloat, ss[3].parseFloat)
    mkInt
  template mkProc5(s:untyped):IntProc =
    proc mkInt(T,V:Integrator; steps:int):Integrator {.gensym.} =
      `mk s`(T = T, V = V, steps = steps, ss[1].parseFloat, ss[2].parseFloat, ss[3].parseFloat, ss[4].parseFloat)
    mkInt
  let ss = s.split(',')
  # Omelyan's triple star integrators, see Omelyan et. al. (2003)
  case ss[0]:
  of "2MN":
    if ss.len == 1: return mkProc1(Omelyan2MN)
    else: return mkProc2(Omelyan2MN)
  of "4MN5FP":
    if ss.len == 1: return mkProc1(Omelyan4MN5FP)
    elif ss.len == 2: return mkProc2(Omelyan4MN5FP)
    elif ss.len == 3: return mkProc3(Omelyan4MN5FP)
    elif ss.len == 4: return mkProc4(Omelyan4MN5FP)
    elif ss.len == 5: return mkProc5(Omelyan4MN5FP)
    else: return mkProc2(Omelyan4MN5FP)
  of "4MN5FV":
    if ss.len == 1: return mkProc1(Omelyan4MN5FV)
    elif ss.len == 2: return mkProc2(Omelyan4MN5FV)
    elif ss.len == 3: return mkProc3(Omelyan4MN5FV)
    elif ss.len == 4: return mkProc4(Omelyan4MN5FV)
    elif ss.len == 5: return mkProc5(Omelyan4MN5FV)
    else: return mkProc2(Omelyan4MN5FV)
  of "6MN7FV": return mkProc1(Omelyan6MN7FV)
  of "4MN3F1GP":  # lambda = 0.2725431326761773  is  FUEL f3g a0=0.109
    if ss.len == 1: return mkProc1(Omelyan4MN3F1GP)
    else: return mkProc2(Omelyan4MN3F1GP)
  of "4MN4F2GVG": return mkProc1(Omelyan4MN4F2GVG)
  of "4MN4F2GV": return mkProc1(Omelyan4MN4F2GV)
  of "4MN5F1GV": return mkProc1(Omelyan4MN5F1GV)
  of "4MN5F1GP": return mkProc1(Omelyan4MN5F1GP)
  of "4MN5F2GV": return mkProc1(Omelyan4MN5F2GV)
  of "4MN5F2GP": return mkProc1(Omelyan4MN5F2GP)
  of "6MN5F3GP": return mkProc1(Omelyan6MN5F3GP)
  else:
    qexError "Cannot parse integrator: '", s, "'\n",
      """Available integrators (with default parameters):
      2MN,0.1931833275037836
      4MN5FP,0.2750081212332419,-0.1347950099106792,-0.08442961950707149,0.3549000571574260
      4MN5FV,0.2539785108410595,-0.03230286765269967,0.08398315262876693,0.6822365335719091
      6MN7FV
      4MN3F1GP,0.2470939580390842
      4MN4F2GVG
      4MN4F2GV
      4MN5F1GV
      4MN5F1GP
      4MN5F2GV
      4MN5F2GP
      6MN5F3GP"""

qexinit()

letParam:
  lat = @[8,8,8,8]
  beta = 6.0
  adjFac = -0.25
  tau = 2.0
  trajs = 10
  seed:uint64 = int(1000*epochTime())
  gintalg:IntProc = "4MN5F2GP"
  gsteps = 4
  mass = @[0.1]  # mass for each staggered species
  hmasses0 = @[0.2,0.4]  # Hasenbusch masses for mass[0]
  hmasses1 = if mass.len>1: hmasses0 else: @[]  # Hasenbusch masses for mass[1]
  hmasses2 = if mass.len>2: hmasses0 else: @[]  # Hasenbusch masses for mass[2]
  hmasses3 = if mass.len>3: hmasses0 else: @[]  # Hasenbusch masses for mass[3]
  hmasses4 = if mass.len>4: hmasses0 else: @[]  # Hasenbusch masses for mass[4]
  fintalg:IntProc = "4MN5F2GP"
  fsteps = repeat(4, mass.len)  # nsteps for each mass
  hfsteps0 = fsteps[0].repeat hmasses0.len  # nsteps for Hasenbusch masses 0
  hfsteps1 = fsteps[0].repeat hmasses1.len  # nsteps for Hasenbusch masses 1
  hfsteps2 = fsteps[0].repeat hmasses2.len  # nsteps for Hasenbusch masses 2
  hfsteps3 = fsteps[0].repeat hmasses3.len  # nsteps for Hasenbusch masses 3
  hfsteps4 = fsteps[0].repeat hmasses4.len  # nsteps for Hasenbusch masses 4
  arsq = 1e-20  # CG r^2 for fermion action
  frsq = repeat(1e-12, mass.len)  # CG r^2 for fermion force for each mass
  hfrsq0 = frsq[0].repeat hmasses0.len  # frsq for Hasenbusch masses 0
  hfrsq1 = frsq[0].repeat hmasses1.len  # frsq for Hasenbusch masses 1
  hfrsq2 = frsq[0].repeat hmasses2.len  # frsq for Hasenbusch masses 2
  hfrsq3 = frsq[0].repeat hmasses3.len  # frsq for Hasenbusch masses 3
  hfrsq4 = frsq[0].repeat hmasses4.len  # frsq for Hasenbusch masses 4
  pbpmass = mass
  pbpreps = repeat(1, pbpmass.len)
  pbprsq = arsq
  maxits = 10000
  useFG2:bool = 0

echoParams()

if mass.len > 5:
  qexError "Unimlemented for mass: ", mass
if mass.len != fsteps.len or
    mass.len != frsq.len:
  qexError "Parameters for staggered species mismatch."

let
  hmasses = @[hmasses0, hmasses1, hmasses2, hmasses3, hmasses4][0..<mass.len]
  hfsteps = @[hfsteps0, hfsteps1, hfsteps2, hfsteps3, hfsteps4][0..<mass.len]
  hfrsq = @[hfrsq0, hfrsq1, hfrsq2, hfrsq3, hfrsq4][0..<mass.len]

for k in 0..<mass.len:
  if hmasses[k].len != hfsteps[k].len or
      hmasses[k].len != hfrsq[k].len:
    qexError "Hasenbusch parameters lengths mismatch."

if pbpmass.len != pbpreps.len:
  qexError "The lengths of pbpmass and pbpreps differ."

let
  lo = lat.newLayout
  #gc = GaugeActionCoeffs(plaq:6)
  gc = GaugeActionCoeffs(plaq: beta, adjplaq: beta*adjFac)
  vol = lo.physVol

var r = lo.newRNGField(RngMilc6, seed)
var R:RngMilc6  # global RNG
R.seed(seed, 987654321)

var g = lo.newgauge
#g.random r
g.unit

echo 6.0*g.plaq
echo gc.actionA g

var
  p = lo.newgauge
  f = lo.newgauge
  g0 = lo.newgauge
  sg0 = lo.newgauge

var ftmp = lo.ColorVector()
type CV = typeof(ftmp)
var
  phi = newseq[seq[CV]](mass.len)
  psi = newseq[seq[CV]](mass.len)
for k in 0..<mass.len:
  phi[k] = newseq[typeof(ftmp)](hmasses[k].len+1)
  psi[k] = newseq[typeof(ftmp)](hmasses[k].len+1)
  for i in 0..<phi[k].len:
    phi[k][i] = lo.ColorVector()
    psi[k][i] = lo.ColorVector()

var pbpsp = initSolverParams()
pbpsp.r2req = pbprsq
pbpsp.maxits = maxits

var spa = initSolverParams()
#spa.subsetName = "even"
spa.r2req = arsq
spa.maxits = maxits

var spf = newseq[typeof(spa)](mass.len)
for k in 0..<mass.len:
  spf[k] = initSolverParams()
  #spf[k].subsetName = "even"
  spf[k].r2req = frsq[k]
  spf[k].maxits = maxits
  spf[k].verbosity = 0

var spfh = newseq[typeof(spf)](mass.len)
for k in 0..<mass.len:
  spfh[k] = newseq[typeof(spa)](hfrsq[k].len)
  for i in 0..<spfh[k].len:
    spfh[k][i] = initSolverParams()
    spfh[k][i].r2req = hfrsq[k][i]
    spfh[k][i].maxits = maxits
    spfh[k][i].verbosity = 0

var
  info: PerfInfo
  coef = HypCoefs(alpha1:0.4, alpha2:0.5, alpha3:0.5)
echo "smear = ",coef

var
  sg = lo.newGauge  # smeared gauge
  gg = lo.newgauge  # FG backup gauge
  sgg = lo.newgauge  # FG backup smeared gauge
let
  stag = newStag(sg)
  stag0 = newStag(sg0)
  stagg = newStag(sgg)

proc pbp(stag:any) =
  var ftmp2 = lo.ColorVector()
  for k in 0..<pbpmass.len:
    for i in 0..<pbpreps[k]:
      threads:
        ftmp.u1 r
      stag.solve(ftmp2, ftmp, pbpmass[k], pbpsp)
      threads:
        var pbp = ftmp2.norm2
        threadMaster:
          echo "MEASpbp mass ",pbpmass[k]," : ",0.5*pbpmass[k]*pbp/vol.float

proc fgsave =
  threads:
    for mu in 0..<g.len:
      gg[mu] := g[mu]
proc fgload =
  threads:
    for mu in 0..<g.len:
      g[mu] := gg[mu]

proc smearRephase(g: any, sg: any):auto {.discardable.} =
  tic()
  let smearedForce = coef.smear(g, sg, info)
  toc("smear")
  threads:
    sg.setBC
    threadBarrier()
    sg.stagPhase
  toc("BC & Phase")
  smearedForce

template pnorm2(p2:float) =
  threads:
    var p2t = 0.0
{'fraction_non_alphanumeric': 0.1358678519639367, 'fraction_numerical': 0.06966773847802786, 'mean_word_length': 3.6137289121582317, 'pattern_counts': {'":': 24, '<': 29, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '787185', 'n_tokens_mistral': 7521, 'n_tokens_neox': 6923, 'n_words': 1496}