From 7c715ec3a2bd21fe60c65d9ae98871640234ec24 Mon Sep 17 00:00:00 2001 From: Hertz4 <88074238+Hertz4@users.noreply.github.com> Date: Wed, 26 Jun 2024 18:35:41 -0700 Subject: [PATCH 1/8] Update README.md --- README.md | 44 ++++++++++++++------------------------------ 1 file changed, 14 insertions(+), 30 deletions(-) diff --git a/README.md b/README.md index b999e27..13d2f80 100644 --- a/README.md +++ b/README.md @@ -31,19 +31,18 @@ Z = np.linspace(-25.,25.,26)*np.pi/beta #Matsubara frequencies Delta = 1.0/(1j*Z-0.5) + 2.0/(1j*Z+0.2) + 0.5/(1j*Z+0.7) # Matsubara functions on these frequencies ``` -With `Delta` and `Z`, one first initialize the Matsubara object: +### Hybridization Fitting +The hybridization fitting is handled by the `hybfit` function. ```python -Imfreq_obj = Matsubara(Delta, Z) +from aaadapol import hybfit ``` - -### Hybridization Fitting -There are two choices for doing hybridization fitting. One can either fit with desired accuracy `eps`: +There are two choices for doing hybridization fitting. One can either fit with desired accuracy tolerance `tol`: ```python -bath_energy, bath_hyb = Imfreq_obj.fitting(tol = 1e-6, flag = "hybfit") +bath_energy, bath_hyb, final_error, func = hybfit(Delta, Z, tol=tol) ``` Or fit with specified number of interpolation points `Np`: ```python -bath_energy, bath_hyb = Imfreq_obj.fitting(Np = 4, flag = "hybfit") +bath_energy, bath_hyb, final_error, func = hybfit(Delta, Z, Np = Np) ``` Here `bath_energy` $E$ and `bath_hyb` $V$ are desired quantities of hybridization orbitals. They satisfy @@ -51,36 +50,21 @@ Here `bath_energy` $E$ and `bath_hyb` $V$ are desired quantities of hybridizatio \Delta(\mathrm i \omega_k)_{mn} \approx \sum_l \frac{V_{lm} V_{ln}^*}{\mathrm i\omega_k - E_l}. ``` -In more sophisticated applications, one might need to specify other flags, such as `maxiter`, `cleanflag` and `disp`. See comments in `matsubara.py` for details. +In more sophisticated applications, one might need to specify other flags, such as `maxiter`, `cleanflag` and `disp`. See the documentation for details. One can look at the final error of the hybridization fitting: ```python -print(Imfreq_obj.final_error) +print(final_error) ``` +### Triqs interface -### Analytic continuation - -Similarly, there are two choices for analytic continuation: - +For triqs users, if the Green's function data `delta_triqs` is stored as a `triqs.gf.Gf` object or `triqs.gf.BlockGf` object, then the hybridization fitting could be done using the `hybfit_triqs` function: ```python -greens_function = Imfreq_obj.fitting(tol = 1e-6, flag = "anacont") +from aaadapol import hybfit_triqs +bathhyb, bathenergy, delta_fit, final_error = hybfit_triqs(delta_triqs, tol=tol, maxiter=500, debug=True) ``` -or - -```python -greens_function = Imfreq_obj.fitting(Np = 4, flag = "anacont") -``` - -The output now is a function evaluator for the Green's functions. For example, if one wish to evaluate the Green's function on `wgrid`, one can do: +### Analytic continuation -```python -wgrid = np.linspace(-1,1,100)+1j*0.01 -G_w = greens_function(wgrid) -``` -and then one can plot the spectral function: -```python -import matplotlib.pyplot as plt -plt.plot(wgrid.real, -np.squeeze(G_w).imag/np.pi) -``` +To use this code for analytic continuation is similar, and we refer to the documentation for details. From 8ae2f35586fed2c7a434aa067fbeb192aad2aa04 Mon Sep 17 00:00:00 2001 From: Hertz4 <88074238+Hertz4@users.noreply.github.com> Date: Wed, 26 Jun 2024 18:36:56 -0700 Subject: [PATCH 2/8] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 13d2f80..85007a8 100644 --- a/README.md +++ b/README.md @@ -62,7 +62,7 @@ print(final_error) For triqs users, if the Green's function data `delta_triqs` is stored as a `triqs.gf.Gf` object or `triqs.gf.BlockGf` object, then the hybridization fitting could be done using the `hybfit_triqs` function: ```python from aaadapol import hybfit_triqs -bathhyb, bathenergy, delta_fit, final_error = hybfit_triqs(delta_triqs, tol=tol, maxiter=500, debug=True) +bathhyb, bathenergy, delta_fit, final_error = hybfit_triqs(delta_triqs, tol=tol, debug=True) ``` ### Analytic continuation From 6b3b25c3faa3732e0fb279fa3d0e6fc42151a9d3 Mon Sep 17 00:00:00 2001 From: Hertz4 <88074238+Hertz4@users.noreply.github.com> Date: Wed, 26 Jun 2024 18:38:09 -0700 Subject: [PATCH 3/8] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 85007a8..51a5f79 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ ## Introduction [`AAAdapol`](https://github.com/Hertz4/AAAdapol) (pronounced "add a pole") is a python package for fitting Matsubara functions with the following form: ```math -G(\mathrm i \omega_k) = \sum_l \frac{V_{lm} V_{ln}^*}{\mathrm i\omega_k - E_l}. +G(\mathrm i \omega_k) = \sum_l \frac{|v_l\rangle\langle v_l|}{\mathrm i\omega_k - E_l}. ``` AAAdapol is short for **A**ntoulas–**A**nderson **Ad**aptive **pol**e-fitting. From f1107c7e0451dc6886404ae2e1fb23342286da36 Mon Sep 17 00:00:00 2001 From: Hertz4 <88074238+Hertz4@users.noreply.github.com> Date: Wed, 26 Jun 2024 18:38:34 -0700 Subject: [PATCH 4/8] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 51a5f79..e673cb9 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ ## Introduction [`AAAdapol`](https://github.com/Hertz4/AAAdapol) (pronounced "add a pole") is a python package for fitting Matsubara functions with the following form: ```math -G(\mathrm i \omega_k) = \sum_l \frac{|v_l\rangle\langle v_l|}{\mathrm i\omega_k - E_l}. +G(\mathrm i \omega_k) = \sum_l \frac{V_lV_l^{\dagger}}{\mathrm i\omega_k - E_l}. ``` AAAdapol is short for **A**ntoulas–**A**nderson **Ad**aptive **pol**e-fitting. From b1f864c96c756ae9840243c191ee057b66e95d30 Mon Sep 17 00:00:00 2001 From: Hertz4 <88074238+Hertz4@users.noreply.github.com> Date: Wed, 26 Jun 2024 18:41:53 -0700 Subject: [PATCH 5/8] Update README.md --- README.md | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index e673cb9..15f6512 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -## Introduction +# Introduction [`AAAdapol`](https://github.com/Hertz4/AAAdapol) (pronounced "add a pole") is a python package for fitting Matsubara functions with the following form: ```math G(\mathrm i \omega_k) = \sum_l \frac{V_lV_l^{\dagger}}{\mathrm i\omega_k - E_l}. @@ -10,7 +10,7 @@ Current applications include We also provide a [TRIQS](https://triqs.github.io/) interface if the Matsubara functions are stored in `triqs` Green's function container. -## Installation +# Installation `AAAdapol` has `numpy` and `scipy` as its prerequisites. [`cvxpy`](https://www.cvxpy.org/) is also required for hybridization fitting of matrix-valued (instead of scalar-valued) Matsubara functions. To install `AAAdapol`, run @@ -19,7 +19,7 @@ pip install aaadapol ``` -## Examples +# Examples In the `examples` directory, we provide two examples [`discrete.ipynb`](https://github.com/Hertz4/AAAdapol/blob/main/example/discrete.ipynb) and [`semicircle.ipynb`](https://github.com/Hertz4/AAAdapol/blob/main/example/semicircle.ipynb), showcasing how to use `AAAdapol` for both discrete spectrum and continuous spectrum. We also demonstrate how to use our code through the triqs interface. Below is a quick introduction through the following toy example: @@ -68,3 +68,10 @@ bathhyb, bathenergy, delta_fit, final_error = hybfit_triqs(delta_triqs, tol=tol, ### Analytic continuation To use this code for analytic continuation is similar, and we refer to the documentation for details. + +# References +To cite this work, ... +### Other References +1. Huang, Zhen, Emanuel Gull, and Lin Lin. "Robust analytic continuation of Green's functions via projection, pole estimation, and semidefinite relaxation." Physical Review B 107.7 (2023): 075151. +2. Mejuto-Zaera, Carlos, et al. "Efficient hybridization fitting for dynamical mean-field theory via semi-definite relaxation." Physical Review B 101.3 (2020): 035143. +3. Nakatsukasa, Yuji, Olivier Sète, and Lloyd N. Trefethen. "The AAA algorithm for rational approximation." SIAM Journal on Scientific Computing 40.3 (2018): A1494-A1522. From d8dde71aed8fba86782e2843bfd9fec827b58974 Mon Sep 17 00:00:00 2001 From: Hertz4 <88074238+Hertz4@users.noreply.github.com> Date: Thu, 27 Jun 2024 02:10:43 -0700 Subject: [PATCH 6/8] Update comments for generating website --- aaadapol/matsubara.py | 347 ++++++++++++++++++++++-------------------- 1 file changed, 180 insertions(+), 167 deletions(-) diff --git a/aaadapol/matsubara.py b/aaadapol/matsubara.py index c64425b..30a7aec 100644 --- a/aaadapol/matsubara.py +++ b/aaadapol/matsubara.py @@ -1,12 +1,12 @@ import numpy as np from .fit_utils import pole_fitting, eval_with_pole, check_psd - -def anacont( +def hybfit( Delta, iwn_vec, tol=None, Np=None, + svdtol=1e-7, solver="lstsq", maxiter=500, mmin=4, @@ -14,74 +14,53 @@ def anacont( verbose=False, ): """ - The main fitting function for both hybridization fitting. + The function for hybridization fitting. + Examples: - -------- + ---------- - func = anacont(Np = Np) # analytic continuation with Np poles - func = anacont(tol = tol) # analytic continuation with fixed error tolerance tol + - Fitting with :math:`N_p` poles: + :code:`hybfit(Delta, iwn_vec, Np = Np)` - Analytic continuation with improved accuracy: - fitting(tol = tol, flag = flag, cleanflag = False) - fitting(Np = Np, flag = flag, cleanflag = False) + - Fitting with fixed error tolerance tol : + :code:`hybfit(Delta, iwn_vec, tol = tol)` Parameters: - -------- - Delta: np.array (Nw, Norb, Norb) - The input Matsubara function in Matsubara frequency - - iwn_vec: np.array (Nw) - The Matsubara frequency vector, complex-valued - - - tol: Fitting error tolreance, float - If tol is specified, the fitting will be conducted with fixed error tolerance tol. - default: None - - Np: number of poles, integer - If Np is specified, the fitting will be conducted with fixed number of poles. - default: None - Np needs to be an odd integer, and number of supoort points is Np + 1. + ------------ + :code:`Delta`: np.array, :math:`(N_w, N_\mathrm{orb}, N_\mathrm{orb})` + The input hybridization function in Matsubara frequency. + :code:`iwn_vec`: np.array, :math:`(N_w)` + The Matsubara frequency vector, complex-valued - solver: string - The solver that is used for optimization. - choices: "lstsq", "sdp" - default: "lstsq" - - maxiter: int - maximum number of iterations - default: 500 - - mmin, mmax: number of minimum or maximum poles, integer - default: mmin = 4, mmax = 50 - if tol is specified, mmin and mmax will be used as the minimum and maximum number of poles. - if Np is specified, mmin and mmax will not be used. + :code:`svdtol`: float, optional + Truncation threshold for bath orbitals while doing SVD of weight matrices in hybridization fitting + default:1e-7 - verbose: bool - whether to display optimization details - default: False + :code:`tol`, :code:`Np`, :code:`solver`, :code:`maxiter`, :code:`mmin`, :code:`mmax`, :code:`verbose`: + see below in anacont + Returns: + --------- - Returns: - -------- - func: function - Analytic continuation function - func(w) = sum_n weight[n]/(w-pol[n]) + :code:`bathenergy` :math:`E`: np.array, :math:`(N_b)` + Bath energy - fitting_error: float - fitting error + :code:`bathhyb` :math:`V`: np.array, :math:`(N_b,N_{\mathrm{orb}})` + Bath hybridization - pol: np.array (Np) - poles obtained from fitting + :code:`final_error`: float + final fitting error - weight: np.array (Np, Norb, Norb) - weights obtained from fitting + :code:`func`: function + Hybridization function evaluator + :math:`f(z) = \sum_n V_{ni}V_{nj}^*/(z-E_n).` """ + # Check dimensions assert len(iwn_vec.shape) == 1 or len(iwn_vec.shape) == 2 if len(iwn_vec.shape) == 2: @@ -110,7 +89,7 @@ def anacont( if Np is not None: pol, weight, fitting_error = pole_fitting( - Delta, wn_vec, Ns=Np+1, maxiter=maxiter, solver=solver, disp=verbose + Delta, wn_vec, Ns=Np + 1, maxiter=maxiter, solver=solver, disp=verbose ) elif tol is not None: pol, weight, fitting_error = pole_fitting( @@ -123,15 +102,20 @@ def anacont( solver=solver, disp=verbose, ) + + bathenergy, bathhyb, bath_mat = obtain_orbitals(pol, weight, svdtol=svdtol) def func(Z): - return eval_with_pole(pol, Z, weight) - return func, fitting_error, pol, weight + return eval_with_pole(bathenergy, Z, bath_mat) + Delta_reconstruct = func(iwn_vec) + final_error = np.max(np.abs(Delta - Delta_reconstruct)) + return bathenergy, bathhyb, final_error, func -def anacont_triqs( +def hybfit_triqs( Delta_triqs, tol=None, Np=None, + svdtol=1e-7, solver="lstsq", maxiter=500, mmin=4, @@ -140,42 +124,54 @@ def anacont_triqs( debug=False ): """ - The triqs interface for analytical continuation. + The triqs interface for hybridization fitting. The function requires triqs package in python. + Examples: - -------- + ---------- + + - Fitting with :math:`N_p` poles: + :code:`hybfit_triqs(delta_triqs, Np = Np)` + + - Fitting with fixed error tolerance tol : + :code:`hybfit_triqs(delta_triqs, tol = tol)` - func = anacont(Np = Np) # analytic continuation with Np poles - func = anacont(tol = tol) # analytic continuation with fixed error tolerance tol Parameters: - -------- - Delta_triqs: triqs Green's function container + ------------ + :code:`Delta_triqs`: triqs Green's function container The input hybridization function in Matsubara frequency - debug: bool + :code:`debug`: bool return additional outputs for debugging. Default: False - tol, Np, cleanflag, maxiter, mmin, mmax, disp: see above in anacont + :code:`svdtol`: float, optional + Truncation threshold for bath orbitals while doing SVD of weight matrices in hybridization fitting + default: 1e-7 + :code:`tol`, :code:`Np`, :code:`solver`, :code:`maxiter`, :code:`mmin`, :code:`mmax`, :code:`verbose`: + same as in hybfit Returns: - -------- - func: function - Analytic continuation function - func(w) = sum_n weight[n]/(w-pol[n]) + --------- - if debug == True: - fitting_error: float - fitting error + :code:`bathhyb`: np.array :math:`(N_b, N_\mathrm{orb})` + Bath hybridization - pol: np.array (Np) - poles obtained from fitting + :code:`bathenergy`: np.array :math:`(N_b,)` + Bath energy - weight: np.array (Np, Norb, Norb) - weights obtained from fitting + :code:`Delta_fit`: triqs Gf or BlockGf + Discretized hybridization function + The input hybridization function in Matsubara frequency + + if debug is True: + :code:`final_error`: float + final fitting error + :code:`weight`: np.array :math:`(N_p, N_\mathrm{orb}, N_\mathrm{orb})` + weights obtained from fitting """ try: from triqs.gf import Gf, BlockGf, MeshImFreq, MeshDLRImFreq @@ -185,39 +181,42 @@ def anacont_triqs( if isinstance(Delta_triqs, Gf) and isinstance(Delta_triqs.mesh, (MeshImFreq, MeshDLRImFreq)): iwn_vec = np.array([iw.value for iw in Delta_triqs.mesh.values()]) - func, fit_error, pol, weight = anacont(Delta_triqs.data, iwn_vec, tol, Np, solver, maxiter, - mmin, mmax, verbose) - print('optimization finished with fitting error {:.3e}'.format(fit_error)) + eps_opt, V_opt, final_error, func = hybfit(Delta_triqs.data, iwn_vec, tol, Np, + svdtol, solver, maxiter, mmin, mmax, verbose) + print('optimization finished with fitting error {:.3e}'.format(final_error)) + + delta_fit = Gf(mesh=Delta_triqs.mesh, target_shape=Delta_triqs.target_shape) + delta_fit.data[:] = func(iwn_vec) if debug: - return func, fit_error, pol, weight + return V_opt.T.conj(), eps_opt, delta_fit, final_error else: - return func + return V_opt.T.conj(), eps_opt, delta_fit elif isinstance(Delta_triqs, BlockGf) and isinstance(Delta_triqs.mesh, (MeshImFreq, MeshDLRImFreq)): - func_list, error_list, pol_list, weight_list = [], [], [], [] + V_list, eps_list, delta_list, error_list = [], [], [], [] for j, (block, delta_blk) in enumerate(Delta_triqs): - func, fit_error, pol, weight = anacont_triqs(delta_blk, tol, Np, solver, maxiter, mmin, - mmax, verbose) - func_list.append(func) + res = hybfit_triqs(delta_blk, tol, Np, svdtol, solver, maxiter, mmin, mmax, verbose, debug) + V_list.append(res[0]) + eps_list.append(res[1]) + delta_list.append(res[2]) if debug: - error_list.append(fit_error) - pol_list.append(pol) - weight_list.append(weight) + error_list.append(res[3]) if debug: - return func_list, error_list, pol_list, weight_list + return V_list, eps_list, BlockGf(name_list=list(Delta_triqs.indices), block_list=delta_list), error_list else: - return func_list + return V_list, eps_list, BlockGf(name_list=list(Delta_triqs.indices), block_list=delta_list) else: raise RuntimeError("Error: Delta_triqs.mesh must be an instance of MeshImFreq or MeshDLRImFreq.") -def hybfit( + + +def anacont( Delta, iwn_vec, tol=None, Np=None, - svdtol=1e-7, solver="lstsq", maxiter=500, mmin=4, @@ -225,51 +224,77 @@ def hybfit( verbose=False, ): """ - The main fitting function for both hybridization fitting. + The function for analytical continuation. + Examples: - -------- + ---------- - hybfit(Delta, iwn_vec, Np = Np) # hybridization fitting with Np poles - hybfit(Delta, iwn_vec, tol = tol) # hybridization fitting with fixed error tolerance tol + - Analytic continuation with :math:`N_p` poles: + :code:`func = anacont(Np = Np)` + + - Fitting with fixed error tolerance tol: + :code:`func = anacont(tol = tol)` + + - Analytic continuation with improved accuracy: + :code:`fitting(Np = Np, flag = flag, solver = "sdp")` Parameters: - -------- - Delta: np.array (Nw, Norb, Norb) - The input hybridization function in Matsubara frequency + ------------ + :code:`Delta`: np.array, :math:`(N_w, N_\mathrm{orb}, N_\mathrm{orb})` + The input hybridization function in Matsubara frequency. - iwn_vec: np.array (Nw) - The Matsubara frequency vector, complex-valued + :code:`iwn_vec`: np.array, :math:`(N_w,)` + The Matsubara frequency vector, complex-valued - svdtol: float, optional - Truncation threshold for bath orbitals while doing SVD of weight matrices in hybridization fitting - default:1e-7 - tol, Np, cleanflag, maxiter, mmin, mmax, disp: see above in anacont + :code:`tol`: Fitting error tolreance, float + If tol is specified, the fitting will be conducted with fixed error tolerance tol. + default: None + :code:`Np`: number of poles, integer + If Np is specified, the fitting will be conducted with fixed number of poles Np. + default: None + Np needs to be an odd integer, and number of supoort points is Np + 1. + :code:`solver`: string + The solver that is used for optimization. + choices: "lstsq", "sdp" + default: "lstsq" - Returns: - -------- + :code:`maxiter`: int + maximum number of iterations + default: 500 + :code:`mmin`, :code:`mmax`: number of minimum or maximum poles, integer + default: mmin = 4, mmax = 50 + if tol is specified, mmin and mmax will be used as the minimum and maximum number of poles. + if Np is specified, mmin and mmax will not be used. - bathenergy: np.array (Nb) - Bath energy + :code:`verbose`: bool + whether to display optimization details + default: False - bathhyb: np.array (Nb, Norb) - Bath hybridization - final_error: float - final fitting error - func: function - Hybridization function evaluator - func(w) = sum_n bathhyb[n,i]*conj(bathhyb[n,j])/(1j*w-bathenergy[n]) + Returns: + --------- + :code:`func`: function + Analytic continuation function: + :math:`f(z) = \sum_n \mathrm{Weight}[n]/(z-\mathrm{pol}[n]).` + + :code:`fitting_error`: float + fitting error + + :code:`pol`: np.array, :math:`(N_p,)` + poles obtained from fitting + + :code:`weight`: np.array, :math:`(N_p, N_\mathrm{orb}, N_\mathrm{orb})` + weights obtained from fitting """ - # Check dimensions assert len(iwn_vec.shape) == 1 or len(iwn_vec.shape) == 2 if len(iwn_vec.shape) == 2: @@ -298,7 +323,7 @@ def hybfit( if Np is not None: pol, weight, fitting_error = pole_fitting( - Delta, wn_vec, Ns=Np + 1, maxiter=maxiter, solver=solver, disp=verbose + Delta, wn_vec, Ns=Np+1, maxiter=maxiter, solver=solver, disp=verbose ) elif tol is not None: pol, weight, fitting_error = pole_fitting( @@ -311,20 +336,15 @@ def hybfit( solver=solver, disp=verbose, ) - - bathenergy, bathhyb, bath_mat = obtain_orbitals(pol, weight, svdtol=svdtol) def func(Z): - return eval_with_pole(bathenergy, Z, bath_mat) - Delta_reconstruct = func(iwn_vec) - final_error = np.max(np.abs(Delta - Delta_reconstruct)) - return bathenergy, bathhyb, final_error, func + return eval_with_pole(pol, Z, weight) + return func, fitting_error, pol, weight -def hybfit_triqs( +def anacont_triqs( Delta_triqs, tol=None, Np=None, - svdtol=1e-7, solver="lstsq", maxiter=500, mmin=4, @@ -333,44 +353,38 @@ def hybfit_triqs( debug=False ): """ - The triqs interface for hybridization fitting. + The triqs interface for analytical continuation. The function requires triqs package in python. - Examples: - -------- - - hybfit_triqs(delta_triqs, Np = Np) # hybridization fitting with Np poles - hybfit_triqs(delta_triqs, tol = tol) # hybridization fitting with fixed error tolerance tol Parameters: - -------- - Delta_triqs: triqs Green's function container + ------------ + :code:`Delta_triqs`: triqs Green's function container The input hybridization function in Matsubara frequency - debug: bool + :code:`debug`: bool return additional outputs for debugging. Default: False - tol, Np, cleanflag, maxiter, mmin, mmax, disp: see above in hybfit + :code:`tol`, :code:`Np`, :code:`solver`, :code:`maxiter`, :code:`mmin`, :code:`mmax`, :code:`verbose`: + same as in anacont - Returns: - -------- - - bathhyb: np.array (Nb, Norb) - Bath hybridization - - bathenergy: np.array (Nb) - Bath energy - - Delta_fit: triqs Gf or BlockGf - Discretized hybridization function - The input hybridization function in Matsubara frequency + Returns: + --------- + :code:`func`: function + Analytic continuation function: + :math:`f(z) = \sum_n \mathrm{Weight}[n]/(z-\mathrm{pol}[n]).` + if debug == True: - final_error: float - final fitting error + :code:`fitting_error`: float + fitting error + + :code:`pol`: np.array, :math:`(N_p,)` + poles obtained from fitting - weight: np.array (Np, Norb, Norb) + :code:`weight`: np.array, :math:`(N_p, N_\mathrm{orb}, N_\mathrm{orb})` weights obtained from fitting + """ try: from triqs.gf import Gf, BlockGf, MeshImFreq, MeshDLRImFreq @@ -380,35 +394,34 @@ def hybfit_triqs( if isinstance(Delta_triqs, Gf) and isinstance(Delta_triqs.mesh, (MeshImFreq, MeshDLRImFreq)): iwn_vec = np.array([iw.value for iw in Delta_triqs.mesh.values()]) - eps_opt, V_opt, final_error, func = hybfit(Delta_triqs.data, iwn_vec, tol, Np, - svdtol, solver, maxiter, mmin, mmax, verbose) - print('optimization finished with fitting error {:.3e}'.format(final_error)) - - delta_fit = Gf(mesh=Delta_triqs.mesh, target_shape=Delta_triqs.target_shape) - delta_fit.data[:] = func(iwn_vec) + func, fit_error, pol, weight = anacont(Delta_triqs.data, iwn_vec, tol, Np, solver, maxiter, + mmin, mmax, verbose) + print('optimization finished with fitting error {:.3e}'.format(fit_error)) if debug: - return V_opt.T.conj(), eps_opt, delta_fit, final_error + return func, fit_error, pol, weight else: - return V_opt.T.conj(), eps_opt, delta_fit + return func elif isinstance(Delta_triqs, BlockGf) and isinstance(Delta_triqs.mesh, (MeshImFreq, MeshDLRImFreq)): - V_list, eps_list, delta_list, error_list = [], [], [], [] + func_list, error_list, pol_list, weight_list = [], [], [], [] for j, (block, delta_blk) in enumerate(Delta_triqs): - res = hybfit_triqs(delta_blk, tol, Np, svdtol, solver, maxiter, mmin, mmax, verbose, debug) - V_list.append(res[0]) - eps_list.append(res[1]) - delta_list.append(res[2]) + func, fit_error, pol, weight = anacont_triqs(delta_blk, tol, Np, solver, maxiter, mmin, + mmax, verbose) + func_list.append(func) if debug: - error_list.append(res[3]) + error_list.append(fit_error) + pol_list.append(pol) + weight_list.append(weight) if debug: - return V_list, eps_list, BlockGf(name_list=list(Delta_triqs.indices), block_list=delta_list), error_list + return func_list, error_list, pol_list, weight_list else: - return V_list, eps_list, BlockGf(name_list=list(Delta_triqs.indices), block_list=delta_list) + return func_list else: raise RuntimeError("Error: Delta_triqs.mesh must be an instance of MeshImFreq or MeshDLRImFreq.") + def obtain_orbitals(pol, weight, svdtol=1e-7): """ obtaining bath orbitals through svd @@ -434,4 +447,4 @@ def check_weight_psd(weight, atol=1e-6): check whether the weight matrices are positive semidefinite """ - return check_psd(weight, atol=atol) + return check_psd(weight, atol=atol) \ No newline at end of file From 25f1a2d97cf7ab0b8e4a974c5f68ee5ca27588a9 Mon Sep 17 00:00:00 2001 From: jasonkaye Date: Fri, 28 Jun 2024 14:27:13 -0400 Subject: [PATCH 7/8] readme edits (AAAdapol -> adapol) --- README.md | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/README.md b/README.md index 15f6512..94ba2df 100644 --- a/README.md +++ b/README.md @@ -1,9 +1,8 @@ -# Introduction -[`AAAdapol`](https://github.com/Hertz4/AAAdapol) (pronounced "add a pole") is a python package for fitting Matsubara functions with the following form: +# adapol: Adaptive Pole Fitting for Quantum Many-Body Physics +[`adapol`](https://github.com/Hertz4/Adapol) (pronounced "add a pole") is a python package for fitting Matsubara functions with the following form: ```math G(\mathrm i \omega_k) = \sum_l \frac{V_lV_l^{\dagger}}{\mathrm i\omega_k - E_l}. ``` -AAAdapol is short for **A**ntoulas–**A**nderson **Ad**aptive **pol**e-fitting. Current applications include (1) hybridization fitting, (2) analytic continuation. @@ -11,16 +10,16 @@ Current applications include We also provide a [TRIQS](https://triqs.github.io/) interface if the Matsubara functions are stored in `triqs` Green's function container. # Installation -`AAAdapol` has `numpy` and `scipy` as its prerequisites. [`cvxpy`](https://www.cvxpy.org/) is also required for hybridization fitting of matrix-valued (instead of scalar-valued) Matsubara functions. +`adapol` has `numpy` and `scipy` as its prerequisites. [`cvxpy`](https://www.cvxpy.org/) is also required for hybridization fitting of matrix-valued (instead of scalar-valued) Matsubara functions. -To install `AAAdapol`, run +To install `adapol`, run ```terminal -pip install aaadapol +pip install adapol ``` # Examples -In the `examples` directory, we provide two examples [`discrete.ipynb`](https://github.com/Hertz4/AAAdapol/blob/main/example/discrete.ipynb) and [`semicircle.ipynb`](https://github.com/Hertz4/AAAdapol/blob/main/example/semicircle.ipynb), showcasing how to use `AAAdapol` for both discrete spectrum and continuous spectrum. We also demonstrate how to use our code through the triqs interface. +In the `examples` directory, we provide two examples [`discrete.ipynb`](https://github.com/Hertz4/adapol/blob/main/example/discrete.ipynb) and [`semicircle.ipynb`](https://github.com/Hertz4/adapol/blob/main/example/semicircle.ipynb), showcasing how to use `adapol` for both discrete spectrum and continuous spectrum. We also demonstrate how to use our code through the triqs interface. Below is a quick introduction through the following toy example: ### Setup @@ -34,7 +33,7 @@ Delta = 1.0/(1j*Z-0.5) + 2.0/(1j*Z+0.2) + 0.5/(1j*Z+0.7) # Matsubara functions o ### Hybridization Fitting The hybridization fitting is handled by the `hybfit` function. ```python -from aaadapol import hybfit +from adapol import hybfit ``` There are two choices for doing hybridization fitting. One can either fit with desired accuracy tolerance `tol`: ```python @@ -61,7 +60,7 @@ print(final_error) For triqs users, if the Green's function data `delta_triqs` is stored as a `triqs.gf.Gf` object or `triqs.gf.BlockGf` object, then the hybridization fitting could be done using the `hybfit_triqs` function: ```python -from aaadapol import hybfit_triqs +from adapol import hybfit_triqs bathhyb, bathenergy, delta_fit, final_error = hybfit_triqs(delta_triqs, tol=tol, debug=True) ``` @@ -70,8 +69,9 @@ bathhyb, bathenergy, delta_fit, final_error = hybfit_triqs(delta_triqs, tol=tol, To use this code for analytic continuation is similar, and we refer to the documentation for details. # References -To cite this work, ... -### Other References +To cite this work, please include a reference to this GitHub repository, and +cite the following references: + 1. Huang, Zhen, Emanuel Gull, and Lin Lin. "Robust analytic continuation of Green's functions via projection, pole estimation, and semidefinite relaxation." Physical Review B 107.7 (2023): 075151. 2. Mejuto-Zaera, Carlos, et al. "Efficient hybridization fitting for dynamical mean-field theory via semi-definite relaxation." Physical Review B 101.3 (2020): 035143. -3. Nakatsukasa, Yuji, Olivier Sète, and Lloyd N. Trefethen. "The AAA algorithm for rational approximation." SIAM Journal on Scientific Computing 40.3 (2018): A1494-A1522. +3. Nakatsukasa, Yuji, Olivier Sète, and Lloyd N. Trefethen. "The AAA algorithm for rational approximation." SIAM Journal on Scientific Computing 40.3 (2018): A1494-A1522. \ No newline at end of file From de35d92e8fcd26b014cabb0b21fe01b054bf6e29 Mon Sep 17 00:00:00 2001 From: Hertz4 Date: Fri, 28 Jun 2024 15:56:41 -0700 Subject: [PATCH 8/8] test commit --- README.md | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/README.md b/README.md index 94ba2df..d114662 100644 --- a/README.md +++ b/README.md @@ -63,7 +63,6 @@ For triqs users, if the Green's function data `delta_triqs` is stored as a `triq from adapol import hybfit_triqs bathhyb, bathenergy, delta_fit, final_error = hybfit_triqs(delta_triqs, tol=tol, debug=True) ``` - ### Analytic continuation To use this code for analytic continuation is similar, and we refer to the documentation for details. @@ -74,4 +73,4 @@ cite the following references: 1. Huang, Zhen, Emanuel Gull, and Lin Lin. "Robust analytic continuation of Green's functions via projection, pole estimation, and semidefinite relaxation." Physical Review B 107.7 (2023): 075151. 2. Mejuto-Zaera, Carlos, et al. "Efficient hybridization fitting for dynamical mean-field theory via semi-definite relaxation." Physical Review B 101.3 (2020): 035143. -3. Nakatsukasa, Yuji, Olivier Sète, and Lloyd N. Trefethen. "The AAA algorithm for rational approximation." SIAM Journal on Scientific Computing 40.3 (2018): A1494-A1522. \ No newline at end of file +3. Nakatsukasa, Yuji, Olivier Sète, and Lloyd N. Trefethen. "The AAA algorithm for rational approximation." SIAM Journal on Scientific Computing 40.3 (2018): A1494-A1522.