Statistical analysis of ensembles of spectral functions

This example is an illustration of the spectral statistical analysis API. In the script below, we use the same input data for a two-component fermionic Green’s function as in another example, but never construct a final solution out of accumulated particular solutions. Instead, we estimate it by computing spectral averages over the ensemble of the particular solutions, along with error bars (dispersions) and two-point correlations. Use of all three resolution functions is shown.

Accumulate particular solutions and compute statistical characteristics

# Import HDFArchive and some TRIQS modules
from h5 import HDFArchive
from import *
import triqs.utility.mpi as mpi

# Import main SOM class
from som import Som
# Import statistical analysis functions
from som.spectral_stats import spectral_avg, spectral_disp, spectral_corr

energy_window = (-4.0, 4.0)  # Energy window to search the solution in

# Parameters for Som.accumulate()
acc_params = {'energy_window': energy_window}
acc_params['verbosity'] = 2  # Verbosity level
acc_params['l'] = 1000       # Number of particular solutions to accumulate
acc_params['f'] = 100        # Number of global updates
acc_params['t'] = 50         # Number of local updates per global update

# Read G(\tau) from an archive.
g_tau = HDFArchive('input.h5', 'r')['g_tau']

# Set the error bars to a constant (all points of g_tau are equally important)
error_bars = g_tau.copy()[:] = 0.01

# Construct a SOM object
cont = Som(g_tau, error_bars, kind="FermionGf")

# Accumulate particular solutions. This may take some time ...

# Compute statistical characteristics of the ensemble of the accumulated
# particular solutions.

# We use a set of energy intervals centered around points of this mesh
w_mesh = MeshReFreq(*energy_window, 401)

for i in (0, 1):
    # Compute spectral averages using different resolution functions.
    avg_rect = spectral_avg(cont, i, w_mesh, "rectangle")
    avg_lorentz = spectral_avg(cont, i, w_mesh, "lorentzian")
    avg_gauss = spectral_avg(cont, i, w_mesh, "gaussian")

    # Compute spectral dispersions using different resolution functions.
    disp_rect = spectral_disp(cont, i, w_mesh, avg_rect, "rectangle")
    disp_lorentz = spectral_disp(cont, i, w_mesh, avg_lorentz, "lorentzian")
    disp_gauss = spectral_disp(cont, i, w_mesh, avg_gauss, "gaussian")

    # Compute two-point correlators using different resolution functions.
    corr_rect = spectral_corr(cont, i, w_mesh, avg_rect, "rectangle")
    corr_lorentz = spectral_corr(cont, i, w_mesh, avg_lorentz, "lorentzian")
    corr_gauss = spectral_corr(cont, i, w_mesh, avg_gauss, "gaussian")

    # On master node, save results to an archive
    if mpi.is_master_node():
        with HDFArchive("results.h5", 'a') as ar:
            gr = ar[str(i)]
            gr['w_mesh'] = w_mesh
            gr['avg_rect'] = avg_rect
            gr['avg_lorentz'] = avg_lorentz
            gr['avg_gauss'] = avg_gauss
            gr['disp_rect'] = disp_rect
            gr['disp_lorentz'] = disp_lorentz
            gr['disp_gauss'] = disp_gauss
            gr['corr_rect'] = corr_rect
            gr['corr_lorentz'] = corr_lorentz
            gr['corr_gauss'] = corr_gauss

Download input file input.h5.

Plot statistical characteristics

from h5 import HDFArchive
from import *
from matplotlib import pyplot as plt
import numpy as np

# Read data from archive
ar = HDFArchive('results.h5', 'r')

w_mesh = ar["0"]["w_mesh"]
w_points = np.array([float(w) for w in w_mesh])

# Plot averaged spectra with error bars derived from estimated dispersions

fig, axes = plt.subplots(3, 2, figsize=(8, 10))
for i in (0, 1):
    gr = ar[str(i)]

    # Rectangular resolution function
    axes[0, i].errorbar(w_points,
    axes[0, i].set_title(r"$A_%d(\omega)$, Rectangular" % i)

    # Lorentzian resolution function
    axes[1, i].errorbar(w_points,
    axes[1, i].set_title(r"$A_%d(\omega)$, Lorentz" % i)

    # Gaussian resolution function
    axes[2, i].errorbar(w_points,
    axes[2, i].set_title(r"$A_%d(\omega)$, Gaussian" % i)

plt.setp(axes, xlim=(w_points[0], w_points[-1]), xlabel=r"$\omega$")

(Source code, png, hires.png, pdf)

# Plot two-point correlators

fig, axes = plt.subplots(3, 2, figsize=(8, 10))
for i in (0, 1):
    gr = ar[str(i)]

    extent = [w_points[0], w_points[-1], w_points[-1], w_points[0]]

    # Rectangular resolution function
    corr = axes[0, i].imshow(gr["corr_rect"], extent=extent, cmap="hot")
    fig.colorbar(corr, ax=axes[0, i])
    axes[0, i].set_title(r"$\sigma_{\omega\omega'}^{(%d)}$, Rectangular" % i)

    # Lorentzian resolution function
    corr = axes[1, i].imshow(gr["corr_lorentz"], extent=extent, cmap="hot")
    fig.colorbar(corr, ax=axes[1, i])
    axes[1, i].set_title(r"$\sigma_{\omega\omega'}^{(%d)}$, Lorentz" % i)

    # Gaussian resolution function
    corr = axes[2, i].imshow(gr["corr_gauss"], extent=extent, cmap="hot")
    fig.colorbar(corr, ax=axes[2, i])
    axes[2, i].set_title(r"$\sigma_{\omega\omega'}^{(%d)}$, Gaussian" % i)


(png, hires.png, pdf)