Full covariance matrix of input data

This example demonstrates three ways to provide information about uncertainty of input data (here, of an imaginary time fermionic Green’s function \(G(\tau)\)),

As estimated error bars \(\sigma(\tau)\);
As a full covariance matrix \(\hat\Sigma_{\tau\tau'}\) ;
As a full covariance matrix \(\hat\Sigma_{\tau\tau'}\) with all eigenvalues shifted up by a constant \(l^2\), where \(l\) is so called filtering level.

Perform analytic continuation

# Import HDFArchive, some TRIQS modules and NumPy
from h5 import HDFArchive
from triqs.gf import *
import triqs.utility.mpi as mpi
import numpy as np

# Import main SOM class and utility functions
from som import Som, fill_refreq, reconstruct

n_w = 801                    # Number of energy slices for G(\omega)
energy_window = (-4.0, 4.0)  # Energy window to search solutions 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']
n_tau = len(g_tau.mesh)

# Estimated error bars on all slices of 'g_tau'
sigma = 0.01

#
# Construct a Som object using the estimated error bars.
#

error_bars = g_tau.copy()
error_bars.data[:] = sigma

cont_eb = Som(g_tau, error_bars, kind="FermionGf")

#
# Construct a Som object using a full covariance matrix.
#

# Create a GF container for the covariance matrix.
# It is defined on an 'n_tau x n_tau' mesh and has a 1D target shape
# with each target element corresponding to one diagonal element of 'g_tau'.
cov_matrix = Gf(mesh=MeshProduct(g_tau.mesh, g_tau.mesh),
                target_shape=(g_tau.target_shape[0],))

# Diagonal covariance matrix equivalent to using 'error_bars'
cov_matrix.data[:, :, 0] = (sigma ** 2) * np.eye(n_tau)
# Add slight correlations between adjacent \tau-slices
cov_matrix.data[:, :, 0] += np.diag(0.00005 * np.ones(n_tau-1), k=1)
cov_matrix.data[:, :, 0] += np.diag(0.00005 * np.ones(n_tau-1), k=-1)

cont_cm = Som(g_tau, cov_matrix, kind="FermionGf")

#
# Construct a Som object using the same covariance matrix and a finite filtering
# level.
#

# Before using the covariance matrix, SOM will shift its eigenvalues up by fl^2.
fl = 0.01
cont_cmfl = Som(g_tau, cov_matrix, kind="FermionGf", filtering_levels=[fl])

#
# Perform analytic continuation
#

for name, cont in (("error_bars", cont_eb),
                   ("cov_matrix", cont_cm),
                   ("cov_matrix_fl", cont_cmfl)):
    cont.accumulate(**acc_params)
    cont.compute_final_solution(verbosity=1)

    # Recover G(\omega) on an energy mesh
    g_w = GfReFreq(window=energy_window, n_points=n_w, indices=g_tau.indices)
    fill_refreq(g_w, cont)

    # G(\tau) reconstructed from the solution
    g_rec_tau = g_tau.copy()
    reconstruct(g_rec_tau, cont)

    # On master node, save results to an archive
    if mpi.is_master_node():
        with HDFArchive("results.h5", 'a') as ar:
            ar.create_group(name)
            gr = ar[name]
            gr['g_tau'] = g_tau
            gr['g_w'] = g_w
            gr['g_rec_tau'] = g_rec_tau

Download input file input.h5.

Plot input and reconstructed imaginary-time GF’s

from h5 import HDFArchive
from triqs.gf import *
from matplotlib import pyplot as plt
from triqs.plot.mpl_interface import oplot

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

g_tau = ar['error_bars']['g_tau']

# Plot input and reconstructed G(\tau)
oplot(g_tau, mode='R', lw=0.8, label=r"$G(\tau)$")
oplot(ar['error_bars']['g_rec_tau'][0, 0], mode='R', lw=0.8,
      label=r"$G^{rec}(\tau)$, error bars")
oplot(ar['cov_matrix']['g_rec_tau'][0, 0], mode='R', lw=0.8,
      label=r"$G^{rec}(\tau)$, covariance matrix")
oplot(ar['cov_matrix_fl']['g_rec_tau'][0, 0], mode='R', lw=0.8,
      label=r"$G^{rec}(\tau)$, covariance matrix with filtering")

plt.xlim((0, g_tau.mesh.beta))
plt.ylabel(r"$G(\tau)$")
plt.legend(loc="lower center")

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

Plot spectral functions

from h5 import HDFArchive
from triqs.gf import *
from matplotlib import pyplot as plt
from triqs.plot.mpl_interface import oplot

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

# Plot spectral functions
oplot(ar['error_bars']['g_w'][0, 0],
      mode='S', lw=0.8, label=r"Error bars")
oplot(ar['cov_matrix']['g_w'][0, 0],
      mode='S', lw=0.8, label=r"Covariance matrix")
oplot(ar['cov_matrix_fl']['g_w'][0, 0],
      mode='S', lw=0.8, label=r"Covariance matrix with filtering")

plt.ylim((0, 0.4))
plt.ylabel(r"$A(\omega)$")
plt.legend(loc="upper left")

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