parse_bBMN.py

#!/usr/bin/python
import os
import sys
import glob
import math
# ------------------------------------------------------------------
# Parse bosonic BMN output files for a single ensemble,
# shuffling the extracted data into dedicated files for plotting
# Normalize the Polyakov loop data by Nc and the bosonic action by Nc^2-1

# First make sure we're calling this from the right place
if not os.path.isdir('Out'):
  print("ERROR: Out/ does not exist")
  sys.exit(1)

ERRFILE = open('ERRORS', 'w')
MISSINGFILES = open('MISSING', 'w')

# Physical observables
SB = open('data/SB.csv', 'w')
print("MDTU,S_B", file=SB)
MYERS = open('data/Myers.csv', 'w')
print("MDTU,Myers/Nt", file=MYERS)
RATIO = open('data/ratio.csv', 'w')
print("MDTU,Ratio(SO6/SO3)", file=RATIO)
EXTENT = open('data/extent.csv', 'w')
print("MDTU,Extent", file=EXTENT)
ENERGY = open('data/energy.csv', 'w')
print("MDTU,E/N^2,E_prime", file=ENERGY)
POLY = open('data/poly.csv', 'w')
print("ReTr(L),ImTr(L)", file=POLY)
POLY_MOD = open('data/poly_mod.csv', 'w')
print("MDTU,|Tr(L)|,ReTr(L),ImTr(L)", file=POLY_MOD)
PL_EIG = open('data/PLeig.csv', 'w')
SCALAR_SQUARES = open('data/scalarsquares.csv', 'w')
print("MDTU,Tr(X1^2),Tr(X2^2),Tr(X3^2),Tr(X4^2),Tr(X5^2),Tr(X6^2),Tr(X7^2),Tr(X8^2),Tr(X9^2)", file=SCALAR_SQUARES)
SCALAR_EIG_AVE = open('data/scalar_eig_ave.csv', 'w')
print("MDTU,min_ave,max_ave", file=SCALAR_EIG_AVE)
SCALAR_EIG = open('data/scalar_eig.csv', 'w')
print("MDTU,min_min,min_max,max_min,max_max", file=SCALAR_EIG)
SCALAR_EIG_WIDTHS = open('data/scalar_eig_widths.csv', 'w')
print("MDTU,min_width,max_width", file=SCALAR_EIG_WIDTHS)

# Evolution observables
ACCP = open('data/accP.csv', 'w')
print("t,accP", file=ACCP)
EXP_DS = open('data/exp_dS.csv', 'w')
print("t,e^(-dS)", file=EXP_DS)
DELTAS = open('data/deltaS.csv', 'w')
print("t,deltaS", file=DELTAS)
ABS_DS = open('data/abs_dS.csv', 'w')
print("t,|deltaS|", file=ABS_DS)
FORCE = open('data/force.csv', 'w')
print("t,B", file=FORCE)
WALLTIME = open('data/walltime.csv', 'w')
print("t,walltime", file=WALLTIME)
WALLTU = open('data/wallTU.csv', 'w')
print("t,cost", file=WALLTU)

# Run parameters
NSTEP = open('data/Nstep.csv', 'w')
print("t,N_f,N_g", file=NSTEP)
STEPSIZE = open('data/stepsize.csv', 'w')
print("t,eps_f,eps_g", file=STEPSIZE)
TLENGTH = open('data/tlength.csv', 'w')
print("t,L", file=TLENGTH)
KEY = open('data/key.csv', 'w')
print("t,file", file=KEY)
TU = open('data/TU.csv', 'w')
print("t,MDTU", file=TU)

# Status checks and running sums for the ensemble as a whole
bAct = [-1.0, -1.0]
Myers = [-1.0, -1.0]
ratio = [-1.0, -1.0]
energy = [-1.0, -1.0]
prime = [-1.0, -1.0]
oldcfg = 0
oldstamp = "start"
traj = 0
MDTU = 0
# ------------------------------------------------------------------



# ------------------------------------------------------------------
# Cycle through files out.$load-$save from list.txt
for temp_tag in open('list.txt'):
  tag = temp_tag.rstrip()
  # Check to make sure that actual files are present
  if "Configs" in tag:
    print("No Out/out.* files found")
    print("No Out/out.* files found", file=ERRFILE)
    sys.exit(1)
  load, cfg = tag.split('-')

  # Initialize running sums and set dummy walltime
  # If the walltime isn't overwritten, then the run died
  # or its output file is corrupted
  walltime = -1
  stamp = "start"

  # Open file
  # If not found, move on to next file instead of killing whole program,
  # but print error message so I know there is a problem
  infile = 'Out/out.' + tag
  if not os.path.isfile(infile):
    print("Problem opening", infile)
    print("Problem opening", infile, file=ERRFILE)
    continue    # Skip to next file
#  print infile  # Monitor running status

  # If not starting from first file in this ensemble,
  # or if we seem to have skipped a file,
  # guess approximate starting trajectory
  traj_per_file = -1
  Nt = -1
  for line in open(infile):
    if line.startswith('PLACEHOLDER'):
      # Placeholder file -- error has been addressed as well as possible,
      # but don't print nonsense wall clock time
      walltime = -2

    # Extract Nc for bosonic action and Polyakov loop normalizations
    # Define DIMF = Nc**2 - 1
    elif line.startswith('BMN, '):
      temp = line.split(',')
      Nc = float(((temp[1]).split())[2])
      DIMF = Nc**2 - 1.0

      # Set up PL_EIG header if this is the first file
      if traj == 0:
        toprint = ''
        for i in range(int(Nc)):
          toprint += ',eig' + str(i)
        print("MDTU%s" % (toprint), file=PL_EIG)

    # Extract Nt, lambda_lat and mu_lat for other normalizations
    elif line.startswith('nt '):
      Nt = float((line.split())[1])

    elif line.startswith('trajecs '):
      traj_per_file = int((line.split())[1])
      endtraj = traj + traj_per_file

    elif line.startswith('lambda '):
      lambda_lat = float((line.split())[1])
      cube_root = lambda_lat**(1.0 / 3.0)
      norm_extent = 0.5 / (cube_root * cube_root)
      norm_energy = 1.0 / (Nt * Nc * Nc * cube_root)
      norm_prime = 1.0 / (Nt * Nt * cube_root * cube_root)

    elif line.startswith('mu '):
      mu_lat = float((line.split())[1])
      norm_Myers = 1.0 / (mu_lat * Nc * Nc)
      break       # Don't go through whole file yet

  if traj_per_file < 0:
    print(infile, "never defines number of trajectories")
    print(infile, "never defines number of trajectories", file=ERRFILE)
    continue    # Skip to next file
  elif Nt < 0:
    print(infile, "never defines Nt")
    print(infile, "never defines Nt", file=ERRFILE)
    continue    # Skip to next file

  if (traj == 0 and int(load) > 0) or (int(load) != oldcfg):
    print(infile, "misses MDTU before", load)
    traj = int(load)
    endtraj = traj + traj_per_file
  # ----------------------------------------------------------------



  # ----------------------------------------------------------------
  # Cycle through lines in the file
  goodtogo = True
  for line in open(infile):
    # Check for premature termination (e.g., layout problem)
    if line.startswith('termination'):
      print(infile, "reports premature termination")
      print(infile, "reports premature termination", file=ERRFILE)
      goodtogo = False
      break

    # Check for unitarity and anti-hermiticity problems
    # Report for info, but don't skip if found
    if line.startswith('Unitarity problem'):
      print(infile, "reports unitarity problem")
      print(infile, "reports unitarity problem", file=ERRFILE)
      break
    if line.startswith('Anti-hermiticity problem'):
      print(infile, "reports anti-hermiticity problem")
      print(infile, "reports anti-hermiticity problem", file=ERRFILE)
      break

  if not goodtogo:
    continue        # Skip this file, which appears malformed

  # At this point we should be able to begin
  oldcfg = int(cfg)
  SCALAR_SQUARES_READY = -1     # To skip duplicate
  scalar_eig_ave = ''
  scalar_eig_ext = ''
  scalar_eig_width = ''
  for line in open(infile):
    # Extract constant run parameters
    if line.startswith('traj_length '):
      tlength = float((line.split())[1])
      print("%d,%g" % (endtraj, tlength), file=TLENGTH)
    elif line.startswith('nstep '):
      Nstep = int((line.split())[1])
      stepsize = tlength / float(Nstep)
    elif line.startswith('nstep_gauge '):
      Nstep_gauge = int((line.split())[1])
      stepsize_gauge = stepsize / float(2.0 * Nstep_gauge)
      Nstep_gauge *= 2 * Nstep
      print("%d,%d,%d" % (endtraj, Nstep, Nstep_gauge), file=NSTEP)
      print("%d,%g,%g" % (endtraj, stepsize, stepsize_gauge), file=STEPSIZE)

    elif line.startswith('Machine '):
      cpus = int((line.split())[5])

    # Check that the file loaded the appropriate configuration
    elif line.startswith('Time stamp '):
      if stamp == "start":    # Loading configuration
        stamp = line.rstrip()
        if stamp != oldstamp and oldstamp != "start":
          print(infile, "time stamp doesn't match final", oldstamp)
          print(infile, "time stamp doesn't match final", oldstamp, file=ERRFILE)
      else:                   # Saving configuration
        oldstamp = line.rstrip()    # Prepare for next file
        stamp = "start"
    # ------------------------------------------------------------

    # ------------------------------------------------------------
    # Now extract evolution observables and physical observables
    # Acceptance comes before (most) measurements
    # The exception is the action
    # This is always measured twice (before and after the trajectory)
    # and which value we want depends on the accept/reject step...
    # Normalize using Nt and DIMF = Nc**2 - 1
    elif line.startswith('action: so3 '):
      if bAct[0] < 0:
        split = line.split()
        bAct[0] = float(split[10]) / (DIMF * Nt)
        td = [float(split[2]), float(split[4]), \
              float(split[6]), float(split[8])]
        Myers[0] = td[3] / Nt
        ratio[0] = td[1] / td[0]
        energy[0] = 3.0 * td[2] + 2.0 * (td[0] + td[1]) + 2.5 * td[3]
        prime[0] = 6.0 * td[2] + 2.0 * (td[0] + td[1]) + 3.75 * td[3]
      elif bAct[1] < 0:
        split = line.split()
        bAct[1] = float(split[10]) / (DIMF * Nt)
        td = [float(split[2]), float(split[4]), \
              float(split[6]), float(split[8])]
        Myers[1] = td[3] / Nt
        ratio[1] = td[1] / td[0]
        energy[1] = 3.0 * td[2] + 2.0 * (td[0] + td[1]) + 2.5 * td[3]
        prime[1] = 6.0 * td[2] + 2.0 * (td[0] + td[1]) + 3.75 * td[3]
      else:
        print(infile, "lists too many action computations")
        print(infile, "lists too many action computations", file=ERRFILE)

    elif ' delta S = ' in line:
      SCALAR_SQUARES_READY = 1
      traj += 1
      temp = MDTU + tlength
      MDTU = round(temp, 3)   # Round off digits in MDTU
      print("%d,%g" % (traj, MDTU), file=TU)
      print("%g,%s" % (traj, tag), file=KEY)

      dS = (line.split())[4]
      # Strip nasty comma from dS, if necessary
      if ',' in dS:
        dS = (dS.split(','))[0]
      print("%d,%g" % (traj, float(dS)), file=DELTAS)
      print("%d,%g" % (traj, math.exp(-1.0 * float(dS))), file=EXP_DS)

      # For RMS, don't have an easy way to average over many measurements
      # Instead just print out absolute value and consider its running average
      print("%d,%g" % (traj, abs(float(dS))), file=ABS_DS)

      # Acceptance is smeared out by running averages
      if line.startswith('ACCEPT'):
        print("%d,1" % traj, file=ACCP)
        print("%g,%g" % (MDTU, bAct[1]), file=SB)
        Myers[1] *= norm_Myers
        print("%g,%g" % (MDTU, Myers[1]), file=MYERS)
        print("%g,%g" % (MDTU, ratio[1]), file=RATIO)
        energy[1] *= norm_energy
        prime[1] *= norm_prime
        print("%g,%g,%g" % (MDTU, energy[1], prime[1]), file=ENERGY)
      else:
        print("%d,0" % traj, file=ACCP)
        print("%g,%g" % (MDTU, bAct[0]), file=SB)
        Myers[0] *= norm_Myers
        print("%g,%g" % (MDTU, Myers[0]), file=MYERS)
        print("%g,%g" % (MDTU, ratio[0]), file=RATIO)
        energy[0] *= norm_energy
        prime[0] *= norm_prime
        print("%g,%g,%g" % (MDTU, energy[0], prime[0]), file=ENERGY)
      bAct = [-1.0, -1.0]                          # Reset
      Myers = [-1.0, -1.0]
      ratio = [-1.0, -1.0]
      energy = [-1.0, -1.0]
      prime = [-1.0, -1.0]

    # Forces -- take maxima rather than average if possible
    elif line.startswith('MONITOR_FORCE_GAUGE '):
      force_g = float((line.split())[-1])
      print("%d,%g" % (traj, force_g), file=FORCE)
    # ------------------------------------------------------------

    # ------------------------------------------------------------
    # Polyakov loop and PL eigenvalues
    # First normalized using Nc extracted above
    elif line.startswith('GMES '):
      temp = line.split()
      poly_r = float(temp[1]) / Nc
      poly_i = float(temp[2]) / Nc
      print("%g,%g" % (poly_r, poly_i), file=POLY)
      poly_mod = math.sqrt(poly_r**2 + poly_i**2)
      print("%g,%g,%g,%g" % (MDTU, poly_mod, poly_r, poly_i), file=POLY_MOD)

    # Format: LINES_EIG {Nc x phase}
    # Check that all phases are within [-pi, pi),
    # accounting for rounding in output files
    # Can be commented out to speed up analysis
    elif line.startswith('LINES_EIG '):
      PL_EIG.write("%g," % (MDTU))
      temp = line.split()
      for i in range(int(Nc - 1)):
        phase = float(temp[1 + i])
        if phase > 3.142 or phase < -3.142:
          print(infile, "phase %.4g exceeds [-pi, pi)" % (phase))
          print(infile, "phase %.4g exceeds [-pi, pi)" % (phase), file=ERRFILE)
        PL_EIG.write("%g," % (phase))
      print("%g" % (float(temp[-1])), file=PL_EIG)
    # ------------------------------------------------------------

    # ------------------------------------------------------------
    # Nine scalar squares
    # Skip duplicate measurement before first trajectory
    elif line.startswith('SCALAR SQUARES ') and SCALAR_SQUARES_READY > 0:
      temp = line.split()
      X1 = float(temp[2])
      X2 = float(temp[3])
      X3 = float(temp[4])
      X4 = float(temp[5])
      X5 = float(temp[6])
      X6 = float(temp[7])
      X7 = float(temp[8])
      X8 = float(temp[9])
      X9 = float(temp[10])
      print("%g,%g,%g,%g,%g,%g,%g,%g,%g,%g" \
                               % (MDTU, X1, X2, X3, X4, X5, X6, X7, X8, X9), file=SCALAR_SQUARES)
      extent = norm_extent * (X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + X9)
      print("%g,%g" % (MDTU, extent), file=EXTENT)
    # ------------------------------------------------------------

    # ------------------------------------------------------------
    # Scalar eigenvalues
    # Only look at largest and smallest (most negative) to keep this manageable
    elif line.startswith('SCALAR_EIG '):
      temp = line.split()
      index = int(temp[1])
      if index == 0 or index == Nc - 1:
        scalar_eig_ave += ',' + str(temp[2])
        scalar_eig_ext += ',' + str(temp[4]) + ',' + str(temp[5])
        scalar_eig_width += ',' + str(temp[3])
      if index == Nc - 1:
        print("%g%s" % (MDTU, scalar_eig_ave), file=SCALAR_EIG_AVE)
        print("%g%s" % (MDTU, scalar_eig_ext), file=SCALAR_EIG)
        print("%g%s" % (MDTU, scalar_eig_width), file=SCALAR_EIG_WIDTHS)
        scalar_eig_ave = ''
        scalar_eig_ext = ''
        scalar_eig_width = ''
    # ------------------------------------------------------------

    # Store total walltime to average at the end
    elif line.startswith('Time = '):
      walltime = float((line.split())[2])
  # ----------------------------------------------------------------



  # ----------------------------------------------------------------
  # Check to see if run seems to have finished properly
  if walltime == -1:
    print(infile, "didn't print final timing")
    print(infile, "didn't print final timing", file=ERRFILE)
  elif walltime == -2:
    # Placeholder file -- error has been addressed as well as possible,
    # but don't print nonsense wall clock time
    bAct = [-1.0, -1.0]                          # Reset
    pass
  else:   # We are good to go
    ave_time = walltime / traj_per_file
    print("%d,%g" % (traj, ave_time), file=WALLTIME)
    TUtime = walltime * float(cpus) / (60.0 * tlength * traj_per_file)
    print("%d,%g" % (traj, TUtime), file=WALLTU)
  # ----------------------------------------------------------------
# ------------------------------------------------------------------



# ------------------------------------------------------------------
# Clean up and close down
ERRFILE.close()
MISSINGFILES.close()
SB.close()
MYERS.close()
RATIO.close()
EXTENT.close()
ENERGY.close()
POLY.close()
POLY_MOD.close()
PL_EIG.close()
SCALAR_SQUARES.close()
SCALAR_EIG_AVE.close()
SCALAR_EIG.close()
SCALAR_EIG_WIDTHS.close()
ACCP.close()
EXP_DS.close()
DELTAS.close()
ABS_DS.close()
FORCE.close()
WALLTIME.close()
WALLTU.close()
NSTEP.close()
STEPSIZE.close()
TLENGTH.close()
KEY.close()
TU.close()
# ------------------------------------------------------------------