diff --git a/final_plots.py b/final_plots.py
index 0150b1dedff75573fcf38e3c3ea7843a5a75015b..8b05af296f656ae21b7171824254232cc7817fab 100644
--- a/final_plots.py
+++ b/final_plots.py
@@ -15,6 +15,7 @@ import argparse
import numpy as np
from scipy.interpolate import bisplrep, bisplev, splrep, BSpline, griddata, interp1d
from scipy.special import jn
+from imageio import imwrite
import settings
from data_management import DataManager, KondoImport
from kondo import gen_init_matrix
@@ -535,6 +536,41 @@ def prepare_plotly_csv():
reduced.to_csv("html/vdc0.csv", columns=reduction_dict.values())
+def prepare_RGconvergence():
+ dc_indices = np.array([0, 8, 20, 30, 90, 165])
+ ac_indices = np.array([30, 50, 80, 125, 170])
+
+ real_cmap = lambda x: plt.cm.seismic((x+1)/2)[...,:3]
+ imwrite(f"figdata/colorbar_seismic.png", np.array(0xffff*plt.cm.seismic(np.linspace(0, 1, 0x1000))[...,:3], dtype=np.uint16).reshape((1,0x1000,3)), format="PNG-FI")
+ def export_img(name, array):
+ imwrite(f"figdata/{name}.png", np.array(0xffff*real_cmap(array[::-1]), dtype=np.uint16), format="PNG-FI")
+
+ data = np.load("figdata/omega5_interp.npz")
+ data_vb7 = np.load("figdata/omega5_interp_vb7.npz")
+ # O2
+ omega = 16.5372
+ omega_o2 = 48.695054
+ vdc = data["vdc"][0] / omega
+ vac = data["vac"][:,0] / omega
+
+ for method in ("o2", "o3", "o3a"):
+ gdiff = np.pi*(data[f"gdc_mu_{method}"] - data["gdc_mu_o3p"])
+ normalize = np.abs(gdiff).max()
+ export_img(f"convergence_{method}_diff", gdiff/normalize)
+ print(f"Color scale for {method} diff: {normalize}")
+
+ for method in ("o2", "o3", "o3a", "o3p"):
+ g = np.pi*data[f"gdc_mu_{method}"]
+ with open(f"figdata/convergence_mu_{method}_dc.dat", "w") as file:
+ file.write("vac " + " ".join("gdc%.2f"%vdc[i] for i in dc_indices) + "\n")
+ np.savetxt(file, np.array([vac, *(g[:,i] for i in dc_indices)]).T)
+ with open(f"figdata/convergence_mu_{method}_ac.dat", "w") as file:
+ file.write("vdc " + " ".join("gdc%.2f"%vac[i] for i in ac_indices) + "\n")
+ np.savetxt(file, np.array([vdc, *(g[i] for i in ac_indices)]).T)
+ print("DC voltage: " + " ".join("%.6g"%vdc[i] for i in dc_indices))
+ print("AC voltage: " + " ".join("%.6g"%vac[i] for i in ac_indices))
+
+
if __name__ == "__main__":
from sys import argv
globals()[argv[1]]()
diff --git a/gen_data.py b/gen_data.py
index 1ef10e027892421a51a4871de94a10c6e924902e..dd1acae02dd6fb156c0669980b5e331a31600f8c 100644
--- a/gen_data.py
+++ b/gen_data.py
@@ -18,7 +18,7 @@ from data_management import DataManager
from kondo import Kondo
-def get_ref_nmax(dm, omega, vdc, vac, preset="precise"):
+def get_ref_nmax(dm, omega, vdc, vac, preset="precise", dc_shift=0, ac_shift=0):
parameters = dict(
old = dict(omega=omega, solver_tol_rel=1e-8, solver_tol_abs=1e-10, d=1e9, method="mu", padding=0, good_flags=0x0000, bad_flags=0x1ffc, xL=0.5, voltage_branches=4),
precise = dict(omega=omega, solver_tol_rel=1e-8, solver_tol_abs=1e-10, d=1e9, method="mu", padding=0, good_flags=0x1000, bad_flags=0x0ffc, xL=0.5, voltage_branches=4),
@@ -28,7 +28,7 @@ def get_ref_nmax(dm, omega, vdc, vac, preset="precise"):
assert data.size > 0
vdc_arr, vac_arr = np.broadcast_arrays(vdc, vac)
nmax = -np.ones(vdc_arr.shape, dtype=np.int64).reshape((-1,))
- for i, (vdc, vac) in enumerate(zip(vdc_arr.flat, vac_arr.flat)):
+ for i, (vdc, vac) in enumerate(zip((vdc_arr+dc_shift).flat, (vac_arr+ac_shift).flat)):
selection = (np.abs(data.vdc - vdc)<1e-6) & (np.abs(data.vac - vac)<1e-6)
n = data.nmax[selection]
if n.size == 1:
@@ -258,6 +258,10 @@ python -m frtrg_kondo.gen_data \\
method_group.add_argument("--preset", metavar="str", type=str,
choices=("normal", "precise", "old"), default="precise",
help = "Preset for choosing nmax if set explicitly")
+ method_group.add_argument("--preset_dc_shift", metavar="float", type=float, default=0,
+ help = "Shift in Vdc for getting nmax from existing data")
+ method_group.add_argument("--preset_ac_shift", metavar="float", type=float, default=0,
+ help = "Shift in Vac for getting nmax from existing data")
# Numerical parameters concerning Floquet matrices
numerics_group = parser.add_argument_group(title="Numerical parameters")
@@ -349,7 +353,11 @@ python -m frtrg_kondo.gen_data \\
atol = options.pop("atol"),
method = options.pop("solver_method"),
)
- preset = options.pop("preset", "precise")
+ preset_options = dict(
+ preset = options.pop("preset", "precise"),
+ dc_shift = options.pop("preset_dc_shift", 0),
+ ac_shift = options.pop("preset_ac_shift", 0),
+ )
# Detect number of CPUs (if necessary)
if threads == 0:
@@ -365,7 +373,7 @@ python -m frtrg_kondo.gen_data \\
assert np.abs(kondo_options["vdc"] - vdc) < 1e-9
kondo_options["vdc"] = vdc
if kondo_options["nmax"] < 0:
- kondo_options["nmax"] = get_ref_nmax(dm, kondo_options["omega"], kondo_options["vdc"], kondo_options["vac"], preset)
+ kondo_options["nmax"] = get_ref_nmax(dm, kondo_options["omega"], kondo_options["vdc"], kondo_options["vac"], **preset_options)
if kondo_options["nmax"] < 0:
settings.logger.error(f"Invalid value nmax={nmax}: must be ≥0")
continue
@@ -385,7 +393,7 @@ python -m frtrg_kondo.gen_data \\
lock = mp.Lock()
queue = mp.Queue()
# create processes
- processes = [mp.Process(target=save_data_mp, args=(queue, lock, solver_options, filename, include, find_pole, False, preset)) for i in range(threads)]
+ processes = [mp.Process(target=save_data_mp, args=(queue, lock, solver_options, filename, include, find_pole, False, preset_options)) for i in range(threads)]
# start processes
for p in processes:
p.start()
@@ -397,7 +405,7 @@ python -m frtrg_kondo.gen_data \\
queue.put(None)
-def save_data_mp(queue, lock, solver_options, filename, include='all', find_pole=False, overwrite=False, preset="precise"):
+def save_data_mp(queue, lock, solver_options, filename, include='all', find_pole=False, overwrite=False, preset_options={}):
"""
Generate data points in own process and save them to HDF5 file.
Each process owns one DataManager instance.
@@ -413,7 +421,7 @@ def save_data_mp(queue, lock, solver_options, filename, include='all', find_pole
assert np.abs(kondo_options["vdc"] - vdc) < 1e-9
kondo_options["vdc"] = vdc
if kondo_options["nmax"] < 0:
- kondo_options["nmax"] = get_ref_nmax(dm, kondo_options["omega"], kondo_options["vdc"], kondo_options["vac"], preset)
+ kondo_options["nmax"] = get_ref_nmax(dm, kondo_options["omega"], kondo_options["vdc"], kondo_options["vac"], **preset_options)
if kondo_options["nmax"] < 0:
settings.logger.error(f"Invalid value nmax={kondo_options['nmax']}: must be ≥0")
continue
diff --git a/plot.py b/plot.py
index 614054df706a3ff87bd06237f8b0f31feab34f60..012003d53f4e8e39dcd81bbb6fb8584ab7279fdc 100644
--- a/plot.py
+++ b/plot.py
@@ -2301,7 +2301,7 @@ def overview_3d_contours(dm, show_contours=False, show_surfaces=True, save_image
contour = ax.contour(vdc0_m, vdc0_y, vdc0_z/TK_VOLTAGE, zdir="x", levels=levels, offset=-offset, zorder=10, colors="black")
save_contour_coordinates(contour, "figdata/overview_contour_x.dat", "vdc", "omega", "gdc")
if save_images:
- array = np.array(256*real_cmap(vdc0_m.T[::-1,::-1]), dtype=np.ubyte)
+ array = np.array(255*real_cmap(vdc0_m.T[::-1,::-1]), dtype=np.ubyte)
img = Image.fromarray(array)
img.save("figdata/overview_x.png")
@@ -2315,7 +2315,7 @@ def overview_3d_contours(dm, show_contours=False, show_surfaces=True, save_image
contour = ax.contour(omega5_x, omega5_y, omega5_m, zdir="z", levels=levels, offset=omega5/TK_VOLTAGE+offset, colors="black")
save_contour_coordinates(contour, "figdata/overview_contour_z.dat", "vac", "vdc", "gdc")
if save_images:
- array = np.array(256*real_cmap(omega5_m[::-1]), dtype=np.ubyte)
+ array = np.array(255*real_cmap(omega5_m[::-1]), dtype=np.ubyte)
img = Image.fromarray(array)
img.save("figdata/overview_z.png")
@@ -2329,7 +2329,7 @@ def overview_3d_contours(dm, show_contours=False, show_surfaces=True, save_image
contour = ax.contour(vac0_x, vac0_m, vac0_z/TK_VOLTAGE, zdir="y", levels=levels, offset=-offset, zorder=10, colors="black")
save_contour_coordinates(contour, "figdata/overview_contour_y.dat", "omega", "vac", "gdc")
if save_images:
- array = np.array(256*real_cmap(vac0_m.T[::-1]), dtype=np.ubyte)
+ array = np.array(255*real_cmap(vac0_m.T[::-1]), dtype=np.ubyte)
img = Image.fromarray(array)
img.save("figdata/overview_y.png")
@@ -2374,19 +2374,23 @@ def RGeq_comparison_new(
vdc = data["vdc"][0] / omega
vac = data["vac"][:,0] / omega
- skip = 40
- plot_number_dc = data["vdc"].shape[0] // skip + 1
- plot_number_ac = data["vac"].shape[1] // skip + 1
- vac_max = vac[-1]
- vdc_max = vdc[-1]
- for i in range(plot_number_dc):
- color = plt.cm.viridis(data["vac"][i*skip,0]/(omega*vac_max))
- ax_dc.plot(data["vdc"][i*skip]/omega, cmp[i*skip], color=color)
- ax_dc.plot(data["vdc"][i*skip]/omega, ref[i*skip], color=color, linestyle=":")
- for i in range(plot_number_ac):
- color = plt.cm.viridis(data["vdc"][0,i*skip]/(omega*vdc_max))
- ax_ac.plot(cmp[:,i*skip], data["vac"][:,i*skip]/omega, color=color)
- ax_ac.plot(ref[:,i*skip], data["vac"][:,i*skip]/omega, color=color, linestyle=":")
+ dc_indices = np.array([0, 8, 20, 30, 90, 165])
+ ac_indices = np.array([30, 50, 80, 125, 170])
+ separation = 0.12
+
+ vmax = max(vdc[dc_indices[-1]], vac[ac_indices[-1]])
+
+ #skip = 40
+ #plot_number_dc = data["vdc"].shape[0] // skip + 1
+ #plot_number_ac = data["vac"].shape[1] // skip + 1
+ #for i in range(plot_number_dc):
+ # color = plt.cm.viridis(data["vac"][i*skip,0]/(omega*vac_max))
+ # ax_dc.plot(data["vdc"][i*skip]/omega, cmp[i*skip], color=color)
+ # ax_dc.plot(data["vdc"][i*skip]/omega, ref[i*skip], color=color, linestyle=":")
+ #for i in range(plot_number_ac):
+ # color = plt.cm.viridis(data["vdc"][0,i*skip]/(omega*vdc_max))
+ # ax_ac.plot(cmp[:,i*skip], data["vac"][:,i*skip]/omega, color=color)
+ # ax_ac.plot(ref[:,i*skip], data["vac"][:,i*skip]/omega, color=color, linestyle=":")
#ax_dc.plot(data["vdc"][::skip].T/omega, cmp[::skip].T, color="blue")
#ax_dc.plot(data["vdc"][::skip].T/omega, ref[::skip].T, ":", color="red")
@@ -2402,6 +2406,20 @@ def RGeq_comparison_new(
img = ax_img.imshow(reldiff, cmap=ccmap, norm=plt.Normalize(-reldiff_max, reldiff_max), extent=extent, origin="lower")
#img4 = ax2.imshow(absdiff, cmap=ccmap, norm=plt.Normalize(-absdiff_max, absdiff_max), extent=extent2, origin="lower")
fig.colorbar(img, ax=ax_trash, location="right")
+
+ #ax_img.hlines(vac[ac_indices], 0, 1, transform=ax_img.get_yaxis_transform(), color="black", linewidth=0.5)
+ #ax_img.vlines(vdc[dc_indices], 0, 1, transform=ax_img.get_xaxis_transform(), color="black", linewidth=0.5)
+ for n, i in enumerate(ac_indices[::-1]):
+ color = plt.cm.viridis(vac[i]/vmax)
+ ax_dc.plot(data["vdc"][i]/omega, n*separation + cmp[i], color=color)
+ ax_dc.plot(data["vdc"][i]/omega, n*separation + ref[i], color=color, linestyle=":")
+ ax_img.hlines(vac[i], 0, 1, color=color, transform=ax_img.get_yaxis_transform(), linewidth=0.5)
+ for n, i in enumerate(dc_indices[::-1]):
+ color = plt.cm.viridis(vdc[i]/vmax)
+ ax_ac.plot(n*separation + cmp[:,i], data["vac"][:,i]/omega, color=color)
+ ax_ac.plot(n*separation + ref[:,i], data["vac"][:,i]/omega, color=color, linestyle=":")
+ ax_img.vlines(vdc[i], 0, 1, color=color, transform=ax_img.get_xaxis_transform(), linewidth=0.5)
+
ax_img.set_xlabel(r"$V_\mathrm{avg}~(\Omega)$")
ax_img.set_ylabel(r"$V_\mathrm{osc}~(\Omega)$")
ax_ac.set_xlabel(r"$G~(2e^2/h)$")