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)