Tutorial
Here we walk through some examples found in the latqcdtools/examples
directory. We try to showcase the flexibility of things like the
bootstrap routine, some convenience wrappers
for plotting, and some pre-packaged physics analysis code.
Hadron resonance gas calculation
The hadron resonance gas (HRG) model, and its implementation in the AnalysisToolbox,
is described in some detail here. Below is a small code
that highlights some of the features of the AnalysisToolbox. For example
initialize
simulataneously saves all screen output to the log file HRG.log
along with the git commit hash, so you can track down which version of the
AnalysisToolbox you used.
Here is latqcdtools/examples/main_HRG_simple.py
import numpy as np
import latqcdtools.base.logger as logger
from latqcdtools.physics.HRG import HRG
from latqcdtools.base.readWrite import readTable, writeTable
from latqcdtools.base.initialize import initialize, finalize
# Write terminal output to log file. Includes git commit hash.
initialize('HRG.log')
# Pick a temperature range in MeV
T = np.arange(100, 166, 1)
# Read in hadron names, masses, charges, baryon number, strangeness,
# charm, and degeneracy factor. This table is provided with AnalysisToolbox.
QMHRG_table = '../latqcdtools/physics/HRGtables/QM_hadron_list_ext_strange_2020.txt'
hadrons, M, Q, B, S, C, g = readTable(QMHRG_table, usecols=(0,1,2,3,4,5,6),
dtype="U11,f8,i8,i8,i8,i8,i8")
w = np.array([1 if ba==0 else -1 for ba in B])
# Instantiate HRG object.
QMhrg = HRG(M,g,w,B,S,Q,C)
# This computation is vectorized since T is a numpy array.
logger.info('Computing chi2B.')
chi = QMhrg.gen_chi(T, B_order=2, Q_order=0, S_order=0, C_order=0,
muB_div_T=0.3, muQ_div_T=0, muS_div_T=0, muC_div_T=0)
# Output T and chi2B as columns in this table.
writeTable("chi2B.txt", T, chi, header=['T [MeV]','chi2B (QMHRG)'])
finalize()
In this case we work at fixed \(\mu_B/T=0.3\)
and compute \(\chi_2^B\) as function of temperature. What is needed is as much
relevant input knowledge about hadron bound states as is known; this is
collected by readTable
. This is a wrapper for numpy.loadtxt()
, hence you
see you can pass it many of the same keyword arguments. The next line
for each species whether the gas is bosonic or fermionic.
Finally the HRG
class is instantiated as QMhrg
, which besides
generic conserved charge cumlants like the one computed with gen_chi
,
contains many methods for various thermodynamic observables such as the
pressure and entropy.
The results are saved in a table chi2B.txt
.
Ising model
The AnalysisToolbox is equipped with a Lattice
class. This class takes care of
indexing, assumes you have periodic boundary conditions, and has iterators, i.e.
functions that perform some operation on every site of a subset of the lattice.
This allows you to write your own statistical mechanical simulations.
Here is latqcdtools/examples/main_isingModel.py
import numpy as np
from latqcdtools.physics.lattice import Lattice
from latqcdtools.base.initialize import initialize, finalize
from latqcdtools.base.plotting import latexify, plt, plot_dots, set_params, clearPlot
from latqcdtools.base.printErrorBars import get_err_str
from latqcdtools.statistics.statistics import std_mean
from latqcdtools.base.speedify import parallel_function_eval, DEFAULTTHREADS
from latqcdtools.base.readWrite import writeTable
from latqcdtools.statistics.jackknife import jackknife
import latqcdtools.base.logger as logger
initialize('example_ising.log')
latexify()
#
# Simulation parameters.
#
Nd = 2 # number of dimensions
Tlow = 1.0 # lowest temperature to sample (kB=1)
Thi = 3.0 # highest temperature to sample
h = 0. # external magnetic field
L = 8 # spatial extension
Nequi = 300 # equilibrate with this many MCMC sweeps
Nmeas = 100 # measure this many MCMC sweeps
Nskip = 5 # separate measurements by this many MCMC sweeps
start = 'hot' # start with all spins up (up), down (down), or random (hot)
Tlist = np.linspace(Tlow,Thi,2*DEFAULTTHREADS)
logger.info()
logger.info('Starting parameters:')
logger.info(' Nequi =',Nequi)
logger.info(' Nmeas =',Nmeas)
logger.info(' Nskip =',Nskip)
logger.info(' Latt =',str(L)+'^'+str(Nd))
logger.info(' h =',h)
logger.info()
#
# We wrap our simulation inside a function. This allows us to trivially parallelize our run
# over the temperatures using parallel_function_eval.
#
def runIsingModel(T):
beta = 1/T
# Initialize the random number generator. When carrying out a statistical physic MCMC, it's
# crucially important that you pick a good one. The default_rng() constructor is what
# numpy recommends, which at the time of writing utilizes O'Neill's PCG algorithm.
rng = np.random.default_rng()
# Initialize the lattice object. The first argument says the lattice geometry will be L**Nd,
# while the second argument is an example of the kind of object that will be on each site.
# Here each site has a scalar, so I just put in the number 1.
lat = Lattice( (L,)*Nd, 1 )
def initialize_site(coord):
if start == 'up':
lat.setElement(coord,1)
elif start == 'down':
lat.setElement(coord,-1)
elif start == 'hot':
lat.setElement(coord,rng.choice([1,-1]))
else:
logger.TBError('Unknown start option',start)
# We then initialize all sites. Here interateOverRandom, whose naming convention follows that
# of SIMULATeQCD, is an iterator that applies the function initialize_site to each site.
# The function should take an array-like coord as argument. Look into the Lattice class to
# see what other iterators are possible.
lat.iterateOverRandom(initialize_site)
def getMagnetization():
""" Order parameter is |M|. """
return np.abs(np.sum(lat.grid)/lat.vol)
def mcLocal(coord):
""" The Markov step.
Args:
coord (array-like): local coordinate
"""
s0 = lat.getElement(coord)
NNsum = 0.
for mu in range(Nd):
NNsum += lat.getElement(lat.march(coord,mu,1)) + lat.getElement(lat.march(coord,mu,-1))
dH = 2*beta*s0*NNsum
if dH < 0:
lat.setElement(coord,-s0)
elif rng.random() < np.exp(-dH):
lat.setElement(coord,-s0)
# Equilibrate the lattice by carrying out Nequi equilibration steps.
for iequi in range(Nequi):
lat.iterateOverRandom(mcLocal)
# With the remaining MCMC steps, we measure the magnetization. We skip every Nskip MCMC
# steps to reduce autocorrelation.
magnetizations = []
for imeas in range(Nmeas):
for iskip in range(Nskip):
lat.iterateOverRandom(mcLocal)
magnetizations.append(getMagnetization())
magnetizations = np.array(magnetizations)
def chi_M(M):
""" Magnetic susceptibility. """
return lat.vol*( std_mean(M**2)-std_mean(M)**2 )
def binder(M):
""" Binder cumulant. For the ising model in 2d, this should tend to 2/3 below
Tc, while it should tend to 0 above Tc, in the thermodynamic limit. """
return 1-np.mean(M**4)/(3*np.mean(M**2)**2)
# We have parallelized over the temperatures, so we're not allowed to parallelize over the jackknife:
# that would be nested parallelization!
Mm , Me = jackknife( std_mean, magnetizations, nproc=1 )
chim, chie = jackknife( chi_M , magnetizations, nproc=1 )
Bm , Be = jackknife( binder , magnetizations, nproc=1 )
logger.info('T, <|M|>, chi, B =',round(T,2),get_err_str(Mm,Me),get_err_str(chim,chie),get_err_str(Bm,Be))
return Mm, Me, chim, chie, Bm, Be
# Parallelize over the temperatures. This parallelization strategy is easiest since each temperature
# amounts to an independent run.
data = parallel_function_eval(runIsingModel,Tlist)
res_M = []
res_E = []
res_chi = []
res_chie = []
res_B = []
res_Be = []
for i in range(len(Tlist)):
res_M.append( data[i][0])
res_E.append( data[i][1])
res_chi.append( data[i][2])
res_chie.append(data[i][3])
res_B.append( data[i][4])
res_Be.append( data[i][5])
# Plot the magnetization.
plot_dots(Tlist,res_M,res_E)
set_params(xlabel='$T$',ylabel='$\\ev{|M|}$',title='$V='+str(L)+'^'+str(Nd)+'$ Ising model',alpha_xlabel=0)
plt.show()
clearPlot()
# Plot the magnetic susceptibility.
plot_dots(Tlist,res_chi,res_chie)
set_params(xlabel='$T$',ylabel='$\\ev{\\chi}$',title='$V='+str(L)+'^'+str(Nd)+'$ Ising model',alpha_xlabel=0)
plt.show()
# Record the results in a nicely formatted table.
writeTable('ising_'+str(L)+'_'+str(Nd)+'.d',Tlist,res_M,res_E,res_chi,res_chie,res_B,res_Be,
header=['T','|M|','|M|_err','chi','chi_err','B','Be'])
finalize()
Again, the most useful feature of this example is the Lattice
class. This wraps a
lattice that saves an object of arbitrary type at each site, here scalars. You can move
one site in a direction using march
. In the Metropolis step, we feature the lattice
accessors getElement
and setElement
. To deal with autocorrelation and propagate
error, we use the jackknife class.
Continuum-limit extrapolation
Suppose we want to perform a continuum-limit extrapolation to determine the deconfinement transition temperature \(T_d\) in pure SU(3). The order parameter for this phase transition is given by the Polyakov loop, \(P\). The transition is first-order in the thermodynamic limit, where \(\langle |P|\rangle\) as function of temperature would jump discontinuously at \(T_d\). At finite volume, this abrupt jump becomes smooth, and \(T_d\) is estimated by the inflection point of the curve.
Here is latqcdtools/examples/main_continuumExtrapolate.py
import numpy as np
import latqcdtools.base.logger as logger
from latqcdtools.base.readWrite import readTable
from latqcdtools.base.printErrorBars import get_err_str
from latqcdtools.base.initialize import initialize, finalize
from latqcdtools.math.num_deriv import diff_deriv
from latqcdtools.math.spline import getSpline
from latqcdtools.statistics.statistics import gaudif
from latqcdtools.statistics.bootstr import bootstr_from_gauss
from latqcdtools.physics.continuumExtrap import continuumExtrapolate
from latqcdtools.physics.constants import r0_phys
from latqcdtools.physics.lattice_params import latticeParams
initialize('cont.log')
Nts = [6,8,10,12,14,16,18,20]
Tds, Tderrs = [], []
for Nt in Nts:
T = []
Ns = Nt*3
# Read in Polyakov loop measurements,
beta, PM, PE = readTable('ploop/Nt'+str(Nt)+'.txt',usecols=(0,1,2))
# Create array of temperatures in physical units
for b in beta:
lp = latticeParams(Ns, Nt, b, scaleType='r0')
T.append( lp.getT() )
t = np.linspace(T[0],T[-1],1001)
# Extract Td from inflection point of <|P|> vs T using natural spline
def getTd(pm):
spl = getSpline(T, pm, natural=True)
dPdT = diff_deriv(t, spl)
maxIndex = np.argmax(dPdT)
return t[maxIndex]
# Error in Td estimate comes from 1000 Gaussian bootstrap samples
Td, Tde = bootstr_from_gauss(getTd, PM, PE, 1000)
Tds.append(Td)
Tderrs.append(Tde)
# Perform O(a^4) continuum-limit extrapolation
result, result_err, chidof = continuumExtrapolate( Nts, Tds, Tderrs, order=2, xtype="Nt",
show_results=True, plot_results=True )
# Do a Z-test against literature result,
r0 = r0_phys(year=2014, units="MeVinv")
Tdr0, Tdr0e = r0 * result[0], r0 * result_err[0]
Tdr0_lit, Tdr0_lite = 0.7457, 0.0045
logger.info('q(ours vs. lit) =',gaudif(Tdr0,Tdr0e,Tdr0_lit,Tdr0_lite))
finalize()
Above we show how such an extrapolation is achieved with
the AnalysisToolbox, along with error estimation, plotting the results, and
carrying out a statistical comparison with the known literature value.
We assume you already have results for \(\langle |P|\rangle\) at various \(N_\tau\), which we
read in from tables of the form Nt6.txt
.
For each \(N_\tau\), this code estimates the inflection point of \(\langle |P|\rangle\) as a function
of \(T\) by taking a derivative using a
spline to get \(T_d(N_\tau)\).
Temperatures are calculated in MeV with the help of our
class that collects ensemble parameters. This
procedure is wrapped in a user-defined function getTc
,
so that errors in the \(\langle |P| \rangle\) data can be
conveniently propagated into the error in \(T_d(N_\tau)\)
using a Gaussian bootstrap.
Having the Nts
, Tds
, and Tderrs
, we are ready to perform a
continuum-limit extrapolation. This will perform an
extrapolation to second order in \(a^2\), i.e. \(\mathcal{O}(a^4)\), print the fit
results to screen, and create a plot of the extrapolation for you.
The arrays result
and result_err
contain the best fit parameters
along with their errors, with result[0]
being the continuum value \(T_d\).
In the last block we compare our result with \(T_d r_0\) from the
literature.
The temperatures calculated in this code implicitly had units of
MeV, hence we need \(r_0\) in physical units.
Finally we call gaudif
to carry out a Gaussian difference test or
Z-test, which is implemented in our statistics module.
Statistical bootstrap
Here is latqcdtools/examples/main_bootstrap.py
import numpy as np
from latqcdtools.base.readWrite import readTable
from latqcdtools.math.spline import getSpline
from latqcdtools.base.plotting import plt, plot_dots, plot_band, latexify
from latqcdtools.statistics.bootstr import bootstr_from_gauss
from latqcdtools.base.initialize import initialize, finalize
initialize('example_bootstrap.log')
latexify()
xdata, ydata, edata = readTable("../testing/statistics/wurf.dat", usecols=(0,2,3))
plot_dots(xdata,ydata,edata)
# We will interpolate to these x-values
xspline = np.linspace(np.min(xdata),np.max(xdata),101)
def splineError(data):
""" We assume no errors in the xdata. For the ydata, we will pass the bootstrap
routine ydata along with the errors. The input to this function, data, will then
be generated in each bootstrap bin by drawing normally from ydata with a spread
of edata. In this example we simply get a spline, but you wrap anything you want
inside your bootstrap procedure.
"""
ys = getSpline(xdata,data,3)
return ys(xspline)
ybs, ybserr = bootstr_from_gauss(splineError,data=ydata,data_std_dev=edata,numb_samples=100)
plot_band(xspline,ybs-ybserr,ybs+ybserr,label='bootstrapped interpolation')
plt.show()
finalize()
In the above example we created an interpolated band using the data; there are 100 bootstrap
samples, where each data point is drawn from normal(ydata,edata)
. The error is taken by
default to be the 68-percentile bounds, and the central value is given as the median. Like
the jackknife
from the Ising model example, the bootstrap routine is parallelized by default.