<jupyter_start><jupyter_text>Loading and checking data  Cooling<jupyter_code>data_dir = "../data/TFIM/logsweep/DM/cooling/"

files = sorted(os.listdir(data_dir))

cooling_data = []
for file in files:
    if not file.endswith('.json'): continue
    cooling_data.append(json.load(open(data_dir+file, 'r')))<jupyter_output><empty_output><jupyter_text>density matrix norm check<jupyter_code>print(' L, K, J/B, Log10(Trace of the DM - 1)')
print(
    *sorted((d['L'], 
             d['K'], 
             round(d['J']/d['B'],2), 
             round(np.log10(np.abs(np.sum(d['eigoccs'])-1)), 0)
            ) for d in cooling_data ),
    sep='\n'
)<jupyter_output>L, K, J/B, Log10(Trace of the DM - 1)
(2, 2, 0.2, -5.0)
(2, 2, 1.0, -5.0)
(2, 2, 5.0, -5.0)
(2, 3, 0.2, -4.0)
(2, 3, 1.0, -5.0)
(2, 3, 5.0, -4.0)
(2, 4, 0.2, -4.0)
(2, 4, 1.0, -5.0)
(2, 4, 5.0, -4.0)
(2, 5, 0.2, -5.0)
(2, 5, 1.0, -5.0)
(2, 5, 5.0, -4.0)
(2, 6, 0.2, -4.0)
(2, 6, 1.0, -4.0)
(2, 6, 5.0, -3.0)
(2, 7, 0.2, -5.0)
(2, 7, 1.0, -5.0)
(2, 7, 5.0, -4.0)
(2, 8, 0.2, -4.0)
(2, 8, 1.0, -5.0)
(2, 8, 5.0, -3.0)
(2, 9, 0.2, -3.0)
(2, 9, 1.0, -4.0)
(2, 9, 5.0, -3.0)
(2, 10, 0.2, -6.0)
(2, 10, 1.0, -4.0)
(2, 10, 5.0, -3.0)
(2, 11, 0.2, -4.0)
(2, 11, 1.0, -5.0)
(2, 11, 5.0, -3.0)
(2, 12, 0.2, -3.0)
(2, 12, 1.0, -6.0)
(2, 12, 5.0, -3.0)
(2, 13, 0.2, -6.0)
(2, 13, 1.0, -5.0)
(2, 13, 5.0, -3.0)
(2, 14, 0.2, -3.0)
(2, 14, 1.0, -5.0)
(2, 14, 5.0, -4.0)
(2, 15, 0.2, -3.0)
(2, 15, 1.0, -4.0)
(2, 15, 5.0, -4.0)
(2, 16, 0.2, -4.0)
(2, 16, 1.0, -4.0)
(2, 16, 5.0, -3.0)
(2, 17, 0.2, -3.0)
(2, 17, 1.0, -4.0)
(2, 17, 5.0, -3.0)
(2, 18, 0.2, -3.0)
(2, 18, 1.0, -3.0)
(2, 18, 5.0, -3.0)
(2, 19, 0.2, -3.[...]<jupyter_text>**Note:** the density matrix simulator accumulates numerical errors, producing a non-normalized final density matrix.We cannot get rid of the numerical error, but to get consistent results we nomalize the results (energy, fidelities) during dat analysis.  Reheating<jupyter_code>data_dir = "../data/TFIM/logsweep/DM/reheating/"

files = sorted(os.listdir(data_dir))

reheating_data = []
for file in files:
    if not file.endswith('.json'): continue
    reheating_data.append(json.load(open(data_dir+file, 'r')))<jupyter_output><empty_output><jupyter_text>density matrix norm check<jupyter_code>print(' L, K, J/B, Log10(Trace of the DM - 1)')
print(
    *sorted((d['L'], 
             d['K'], 
             round(d['J']/d['B'],2), 
             round(np.log10(np.abs(np.sum(d['eigoccs'])-1)), 0),
             round(np.sum(d['eigoccs']), 10)
            ) for d in reheating_data ),
    sep='\n'
)<jupyter_output>L, K, J/B, Log10(Trace of the DM - 1)
(2, 2, 0.2, -5.0, 1.0000264645)
(2, 2, 1.0, -5.0, 0.9999960959)
(2, 2, 5.0, -5.0, 0.9999701995)
(2, 3, 0.2, -4.0, 0.9999496937)
(2, 3, 1.0, -5.0, 0.9999920726)
(2, 3, 5.0, -4.0, 0.9998720866)
(2, 4, 0.2, -4.0, 0.9999471009)
(2, 4, 1.0, -6.0, 1.000000596)
(2, 4, 5.0, -4.0, 1.0000449615)
(2, 5, 0.2, -5.0, 1.0000151098)
(2, 5, 1.0, -5.0, 0.9999944866)
(2, 5, 5.0, -4.0, 0.9998138202)
(2, 6, 0.2, -5.0, 1.0000211596)
(2, 6, 1.0, -4.0, 1.0000389218)
(2, 6, 5.0, -3.0, 0.9993733419)
(2, 7, 0.2, -5.0, 0.9999923408)
(2, 7, 1.0, -5.0, 1.0000233501)
(2, 7, 5.0, -4.0, 0.9997420665)
(2, 8, 0.2, -4.0, 1.0001990199)
(2, 8, 1.0, -5.0, 1.0000094771)
(2, 8, 5.0, -3.0, 1.0004398702)
(2, 9, 0.2, -3.0, 0.9994894564)
(2, 9, 1.0, -4.0, 0.9999193847)
(2, 9, 5.0, -3.0, 1.0005251393)
(2, 10, 0.2, -7.0, 1.0000002682)
(2, 10, 1.0, -4.0, 1.0001587272)
(2, 10, 5.0, -3.0, 1.0015248843)
(2, 11, 0.2, -5.0, 0.9999694228)
(2, 11, 1.0, -5.0, 1.0000099838)
(2, 11, 5.0, -3.0, 0.99957622[...]<jupyter_text>Iterative LogSweep cooling<jupyter_code>data_dir = "../data/TFIM/logsweep/DM/iterative/"

files = sorted(os.listdir(data_dir))

iterative_data = []
for file in files:
    if not file.endswith('.json'): continue
    iterative_data.append(json.load(open(data_dir+file, 'r')))<jupyter_output><empty_output><jupyter_text>density matrix norm check<jupyter_code>print(' L, K, J/B, Log10(|Trace of the DM - 1|)')
print(
    *sorted((d['L'], 
             d['K'], 
             round(d['J']/d['B'],2), 
             round(np.log10(np.abs(np.sum(d['eigoccs'])-1)), 0),
             round(np.sum(d['eigoccs']), 10)
            ) for d in iterative_data ),
    sep='\n'
)<jupyter_output>L, K, J/B, Log10(|Trace of the DM - 1|)
(7, 2, 0.2, -3.0, 0.9994182922)
(7, 2, 1.0, -3.0, 0.9991217926)
(7, 2, 5.0, -4.0, 0.9998503257)
(7, 3, 0.2, -3.0, 0.9991117576)
(7, 3, 1.0, -5.0, 0.9999860141)
(7, 3, 5.0, -3.0, 0.9994342903)
(7, 4, 0.2, -3.0, 1.0024900949)
(7, 4, 1.0, -4.0, 1.0002777206)
(7, 4, 5.0, -3.0, 0.9986783145)
(7, 5, 0.2, -3.0, 1.0026191454)
(7, 5, 1.0, -3.0, 0.9977154661)
(7, 5, 5.0, -3.0, 0.999583139)
(7, 6, 0.2, -3.0, 1.0030249655)
(7, 6, 1.0, -2.0, 0.9953251454)
(7, 6, 5.0, -3.0, 0.9985978736)
(7, 7, 0.2, -3.0, 1.0027176132)
(7, 7, 1.0, -2.0, 0.9946759479)
(7, 7, 5.0, -3.0, 0.9978047287)
(7, 8, 0.2, -3.0, 1.0021703085)
(7, 8, 1.0, -2.0, 0.9964862953)
(7, 8, 5.0, -3.0, 0.9977359791)
(7, 9, 0.2, -3.0, 1.0010678535)
(7, 9, 1.0, -2.0, 0.993444512)
(7, 9, 5.0, -3.0, 0.9987197995)
(7, 10, 0.2, -2.0, 0.9943049653)
(7, 10, 1.0, -2.0, 0.9931355926)
(7, 10, 5.0, -3.0, 1.0011528856)
(7, 11, 0.2, -2.0, 0.98604347)
(7, 11, 1.0, -2.0, 0.993367108)
(7, 11, 5.0, -3.0, 0.9979585453[...]<jupyter_text>Available data summary<jupyter_code>print(' K , L, J/B ')
avail_cooling = [(d['K'], d['L'], round(d['J']/d['B'],1)) for d in cooling_data]
avail_reheating = [(d['K'], d['L'], round(d['J']/d['B'],1)) for d in reheating_data]
avail_iterative = [(d['K'], d['L'], round(d['J']/d['B'],1)) for d in iterative_data]

from itertools import product
for K, L, JvB in np.unique(avail_cooling + avail_reheating + avail_iterative, axis=0):
    K = int(K)
    L = int(L)
    if L!=7: continue
    print((K, L, JvB), 
          'C' if (K, L, JvB) in avail_cooling else ' ',
          'R' if (K, L, JvB) in avail_reheating else ' ',
          'It' if (K, L, JvB) in avail_iterative else '  '
         )<jupyter_output>K , L, J/B 
(2, 7, 0.2) C R It
(2, 7, 1.0) C R It
(2, 7, 5.0) C R It
(3, 7, 0.2) C R It
(3, 7, 1.0) C R It
(3, 7, 5.0) C R It
(4, 7, 0.2) C R It
(4, 7, 1.0) C R It
(4, 7, 5.0) C R It
(5, 7, 0.2) C R It
(5, 7, 1.0) C R It
(5, 7, 5.0) C R It
(6, 7, 0.2) C R It
(6, 7, 1.0) C R It
(6, 7, 5.0) C R It
(7, 7, 0.2) C R It
(7, 7, 1.0) C R It
(7, 7, 5.0) C R It
(8, 7, 0.2) C R It
(8, 7, 1.0) C R It
(8, 7, 5.0) C R It
(9, 7, 0.2) C R It
(9, 7, 1.0) C R It
(9, 7, 5.0) C R It
(10, 7, 0.2) C R It
(10, 7, 1.0) C R It
(10, 7, 5.0) C R It
(11, 7, 0.2) C R It
(11, 7, 1.0) C R It
(11, 7, 5.0) C R It
(12, 7, 0.2) C R It
(12, 7, 1.0) C R It
(12, 7, 5.0) C R It
(13, 7, 0.2) C R It
(13, 7, 1.0) C R It
(13, 7, 5.0) C R It
(14, 7, 0.2) C R It
(14, 7, 1.0) C R It
(14, 7, 5.0) C R It
(15, 7, 0.2) C R It
(15, 7, 1.0) C R It
(15, 7, 5.0) C R It
(16, 7, 0.2) C R It
(16, 7, 1.0) C R It
(16, 7, 5.0) C R It
(17, 7, 0.2) C R It
(17, 7, 1.0) C R It
(17, 7, 5.0) C R It
(18, 7, 0.2) C R It
(18, 7, 1.0) C R It
(18, 7, 5.0[...]<jupyter_text>Varying energy gradation number K<jupyter_code>L = 7<jupyter_output><empty_output><jupyter_text>cooling  energy vs K<jupyter_code># L = 7
for JvB in [.2, 1, 5]:
    data_iterator = ((d['K'], d['energy'], np.sum(d['eigoccs'])) 
                     for d in cooling_data 
                     if d['L'] == L and np.isclose(d['J']/d['B'], JvB))
    K_l, E_l, norms_l = zip(*sorted(data_iterator))
    plt.plot(K_l, np.array(E_l)/np.array(norms_l), 'o-', label=f'$J/B={JvB}$')
plt.legend()<jupyter_output><empty_output><jupyter_text>GS infidelity vs K<jupyter_code># L = 7
for JvB in [.2, 1, 5]:
    data_iterator = ((d['K'], gsm_fidelity(d), np.sum(d['eigoccs'])) 
                     for d in cooling_data 
                     if d['L'] == L and np.isclose(d['J']/d['B'], JvB))
    K_l, fidelty_l, norms_l = zip(*sorted(data_iterator))
    fidelty_l /= np.array(norms_l)
    infidelity_l = 1 - np.array(fidelty_l)
    plt.plot(K_l, infidelity_l, 'o-', label=f'$J/B={JvB}$')
plt.legend()
plt.xlabel('K')
plt.ylabel('GS manifold infidelity')
plt.xscale('log')
plt.yscale('log')<jupyter_output><empty_output><jupyter_text>reheating  energy vs K<jupyter_code># L = 7
for JvB in [.2, 1, 5]:
    data_iterator = ((d['K'], d['energy'], np.sum(d['eigoccs'])) 
                     for d in reheating_data 
                     if d['L'] == L and np.isclose(d['J']/d['B'], JvB))
    K_l, E_l, norms_l = zip(*sorted(data_iterator))
    plt.plot(K_l, np.array(E_l)/np.array(norms_l), 'x:', label=f'$J/B={JvB}$')
plt.legend()<jupyter_output><empty_output><jupyter_text>GS infidelity vs K<jupyter_code># L = 7
for JvB in [.2, 1, 5]:
    data_iterator = ((d['K'], gsm_fidelity(d), np.sum(d['eigoccs'])) 
                     for d in reheating_data 
                     if d['L'] == L and np.isclose(d['J']/d['B'], JvB))
    K_l, fidelty_l, norms_l = zip(*sorted(data_iterator))
    fidelty_l /= np.array(norms_l)
    infidelity_l = 1 - np.array(fidelty_l)
    plt.plot(K_l, infidelity_l, 'x:', label=f'$J/B={JvB}$')
plt.legend()
plt.xscale('log')
plt.yscale('log')<jupyter_output><empty_output><jupyter_text>iterative  energy vs K<jupyter_code># L = 7
for JvB in [.2, 1, 5]:
    data_iterator = ((d['K'], d['energy'], np.sum(d['eigoccs'])) 
                     for d in iterative_data 
                     if d['L'] == L and np.isclose(d['J']/d['B'], JvB))
    K_l, E_l, norms_l = zip(*sorted(data_iterator))
    plt.plot(K_l, np.array(E_l)/np.array(norms_l), '+--', label=f'$J/B={JvB}$')
plt.legend()<jupyter_output><empty_output><jupyter_text>GS infidelity vs K<jupyter_code># L = 7
for JvB in [.2, 1, 5]:
    data_iterator = ((d['K'], gsm_fidelity(d), np.sum(d['eigoccs'])) 
                     for d in iterative_data 
                     if d['L'] == L and np.isclose(d['J']/d['B'], JvB))
    K_l, fidelty_l, norms_l = zip(*sorted(data_iterator))
    fidelty_l /= np.array(norms_l)
    infidelity_l = 1 - np.array(fidelty_l)
    plt.plot(K_l, infidelity_l, '+--', label=f'$J/B={JvB}$')
plt.legend()
plt.xscale('log')
plt.yscale('log')<jupyter_output><empty_output><jupyter_text>combined  energy vs K<jupyter_code># L = 7
plt.title(f'TFIM chain $L={L}$. Standard LogSweep density matrix sim')
# cooling
for JvB in [.2, 1, 5]:
    data_iterator = ((d['K'], d['energy'], np.sum(d['eigoccs'])) 
                     for d in cooling_data 
                     if d['L'] == L and np.isclose(d['J']/d['B'], JvB))
    K_l, E_l, norms_l = zip(*sorted(data_iterator))
    plt.plot(K_l, np.array(E_l)/np.array(norms_l), 'o-', label=f'cooling $J/B={JvB}$')
    
plt.gca().set_prop_cycle(None)

# reheating
for JvB in [.2, 1, 5]:
    data_iterator = ((d['K'], d['energy'], np.sum(d['eigoccs'])) 
                     for d in reheating_data 
                     if d['L'] == L and np.isclose(d['J']/d['B'], JvB))
    K_l, E_l, norms_l = zip(*sorted(data_iterator))
    plt.plot(K_l, np.array(E_l)/np.array(norms_l), 'x:', label=f'reheating $J/B={JvB}$')
    
plt.gca().set_prop_cycle(None)

# iterative
for JvB in [.2, 1, 5]:
    data_iterator = ((d['K'], d['energy'], np.sum(d['eigoccs'])) 
                     for d in iterative_data 
                     if d['L'] == L and np.isclose(d['J']/d['B'], JvB))
    K_l, E_l, norms_l = zip(*sorted(data_iterator))
    plt.plot(K_l, np.array(E_l)/np.array(norms_l), '+--', label=f'iterative $J/B={JvB}$')
    
plt.legend()<jupyter_output><empty_output><jupyter_text>logarithmic<jupyter_code># L = 7
plt.figure(figsize=(10, 5))
plt.title(f'TFIM chain $L={L}$. Standard LogSweep density matrix sim')
# cooling
for JvB in [.2, 1, 5]:
    data_iterator = ((d['K'], d['energy'], np.sum(d['eigoccs'])) 
                     for d in cooling_data 
                     if d['L'] == L and np.isclose(d['J']/d['B'], JvB))
    K_l, E_l, norms_l = zip(*sorted(data_iterator))
    plt.plot(K_l, 1 + np.array(E_l)/np.array(norms_l), 'o-', label=f'cooling $J/B={JvB}$')
    
plt.gca().set_prop_cycle(None)

# reheating
for JvB in [.2, 1, 5]:
    data_iterator = ((d['K'], d['energy'], np.sum(d['eigoccs'])) 
                     for d in reheating_data 
                     if d['L'] == L and np.isclose(d['J']/d['B'], JvB))
    K_l, E_l, norms_l = zip(*sorted(data_iterator))
    plt.plot(K_l, 1 + np.array(E_l)/np.array(norms_l), 'x:', label=f'reheating $J/B={JvB}$')
    
plt.gca().set_prop_cycle(None)

# iterative
for JvB in [.2, 1, 5]:
    data_iterator = ((d['K'], d['energy'], np.sum(d['eigoccs'])) 
                     for d in iterative_data 
                     if d['L'] == L and np.isclose(d['J']/d['B'], JvB))
    K_l, E_l, norms_l = zip(*sorted(data_iterator))
    plt.plot(K_l, 1 + np.array(E_l)/np.array(norms_l), '+--', label=f'iterative $J/B={JvB}$')
    
plt.yscale('log')
plt.xscale('log')
plt.ylabel('energy relative to GS '
           r'$(1 + \mathrm{Tr}[\rho H]/\vert E_\mathrm{GS}\vert)$')
plt.legend(bbox_to_anchor = (1, 0.5), loc = 'center left')<jupyter_output><empty_output><jupyter_text>GS infidelity vs K<jupyter_code>L = 7
#plt.title(f'TFIM chain $L={L}$. Standard LogSweep density matrix sim')

JvBlist = [.2, 1, 5]

# cooling
for JvB in JvBlist:
    data_iterator = ((d['K'], gsm_fidelity(d), np.sum(d['eigoccs'])) 
                     for d in cooling_data 
                     if d['L'] == L and np.isclose(d['J']/d['B'], JvB))
    K_l, fidelty_l, norms_l = zip(*sorted(data_iterator))
    fidelty_l /= np.array(norms_l)
    infidelity_l = 1 - np.array(fidelty_l)
    plt.plot(K_l, infidelity_l, '+:', label=f'cooling $J/B={JvB}$')

plt.gca().set_prop_cycle(None)

# reheating
for JvB in JvBlist:
    data_iterator = ((d['K'], gsm_fidelity(d), np.sum(d['eigoccs'])) 
                     for d in reheating_data 
                     if d['L'] == L and np.isclose(d['J']/d['B'], JvB))
    K_l, fidelty_l, norms_l = zip(*sorted(data_iterator))
    fidelty_l /= np.array(norms_l)
    infidelity_l = 1 - np.array(fidelty_l)
    plt.plot(K_l, infidelity_l, 'x:', label=f'reheating $J/B={JvB}$')  

plt.gca().set_prop_cycle(None)

for JvB in JvBlist:
    data_iterator = ((d['K'], gsm_fidelity(d), np.sum(d['eigoccs'])) 
                     for d in iterative_data 
                     if d['L'] == L and np.isclose(d['J']/d['B'], JvB))
    K_l, fidelty_l, norms_l = zip(*sorted(data_iterator))
    fidelty_l /= np.array(norms_l)
    infidelity_l = 1 - np.array(fidelty_l)
    plt.plot(K_l, infidelity_l, 'o:', label=f'iterative $J/B={JvB}$')
    
handles, labels = plt.gca().get_legend_handles_labels()
handles = np.array([[[plt.Line2D([],[],marker='',ls='')] + handles[3*i:3*i+3]] 
                    for i in range(3)]).flatten()
ll = ['$JvB = {JvB}$' for JvB in JvBlist]
labels = ['cooling', *ll, 'reheating', *ll, 'iterative', *ll]
plt.legend(handles, labels, ncol=3, bbox_to_anchor = (0.5, 1), loc='lower center')
plt.yscale('log')
plt.xscale('log')
plt.xlabel('K')
plt.xticks([2, 10, 50], [2, 10, 50])
plt.ylabel('ground space infidelity')
L = 7
#plt.title(f'TFIM chain $L={L}$. Standard LogSweep density matrix sim')

JvBlist = [.2, 1, 5]

def fitfunc(x, a, b):
    return a * x**-b
from scipy.optimize import curve_fit

clrs = plt.rcParams['axes.prop_cycle'].by_key()['color']

# reheating
for JvB, c in zip(JvBlist, clrs):
    data_iterator = ((d['K'], gsm_fidelity(d), np.sum(d['eigoccs'])) 
                     for d in reheating_data 
                     if d['L'] == L and np.isclose(d['J']/d['B'], JvB))
    K_l, fidelty_l, norms_l = zip(*sorted(data_iterator))
    fidelty_l /= np.array(norms_l)
    infidelity_l = 1 - np.array(fidelty_l)
    plt.plot(K_l, infidelity_l, 'x', label=f'reheating $J/B={JvB}$', c=c)  
    
    
# iterative cooling
annotate_xy_list = [
    [4, 0.025],
    [10, 0.18],
    [6, 0.07]
]
for JvB, c, annotate_xy in zip(JvBlist, clrs, annotate_xy_list):
    data_iterator = ((d['K'], gsm_fidelity(d), np.sum(d['eigoccs'])) 
                     for d in iterative_data 
                     if d['L'] == L and np.isclose(d['J']/d['B'], JvB))
    K_l, fidelty_l, norms_l = zip(*sorted(data_iterator))
    fidelty_l /= np.array(norms_l)
    infidelity_l = 1 - np.array(fidelty_l)
    plt.plot(K_l, infidelity_l, '.', label=f'iterative $J/B={JvB}$', color=c)
    
#     popt, pcov = curve_fit(fitfunc, K_l[10:], infidelity_l[10:])
#     plt.plot([2, 50], fitfunc(np.array([2,50]), *popt), '-', c=c, lw=1)
#     plt.annotate(r'$\mathbf{\propto K^{-'+f'{popt[1]:.2f}'+r'}}$', 
#                  annotate_xy, fontsize=16, zorder=10, 
#                  va='center', ha='left')
    
# one-shot cooling
for JvB, c in zip(JvBlist, clrs):
    data_iterator = ((d['K'], gsm_fidelity(d), np.sum(d['eigoccs'])) 
                     for d in cooling_data 
                     if d['L'] == L and np.isclose(d['J']/d['B'], JvB))
    K_l, fidelty_l, norms_l = zip(*sorted(data_iterator))
    fidelty_l /= np.array(norms_l)
    infidelity_l = 1 - np.array(fidelty_l)
    plt.plot(K_l, infidelity_l, '^', label=f'cooling $J/B={JvB}$', c=c)  

    

hndlbl = [[plt.Line2D([],[],marker='',color=clrs[0]), '$J/B=.2$'],
          [plt.Line2D([],[],marker='.',ls='',color='grey'), 'iterative\ncooling'],
          [plt.Line2D([],[],marker='',color=clrs[1]), '$J/B=1$'],
          [plt.Line2D([],[],marker='x',ls='',color='grey'), 'reheating'],
          [plt.Line2D([],[],marker='',color=clrs[2]), '$J/B=5$'],
          [plt.Line2D([],[],marker='^',ls='',color='grey', lw=1), 
           'single-shot\ncooling'],
          ]
plt.legend(*zip(*hndlbl), ncol=3, bbox_to_anchor = (0.5, 1), loc='lower center')
plt.yscale('log')
plt.xscale('log')
plt.xlabel('K')
plt.xticks([2, 10, 50], [2, 10, 50])
plt.ylabel('ground space infidelity')

#plt.savefig('../figures/Ksingleshot.pdf', bbox_inches='tight')<jupyter_output><empty_output><jupyter_text>without cooling<jupyter_code>L = 7
#plt.title(f'TFIM chain $L={L}$. Standard LogSweep density matrix sim')

JvBlist = [.2, 1, 5]

def fitfunc(x, a, b):
    return a * x**-b
from scipy.optimize import curve_fit

clrs = plt.rcParams['axes.prop_cycle'].by_key()['color']

# reheating
for JvB, c in zip(JvBlist, clrs):
    data_iterator = ((d['K'], gsm_fidelity(d), np.sum(d['eigoccs'])) 
                     for d in reheating_data 
                     if d['L'] == L and np.isclose(d['J']/d['B'], JvB))
    K_l, fidelty_l, norms_l = zip(*sorted(data_iterator))
    fidelty_l /= np.array(norms_l)
    infidelity_l = 1 - np.array(fidelty_l)
    plt.plot(K_l, infidelity_l, 'x', label=f'reheating $J/B={JvB}$', c=c)  
    
annotate_xy_list = [
    [4, 0.025],
    [10, 0.18],
    [6, 0.07]
]
for JvB, c, annotate_xy in zip(JvBlist, clrs, annotate_xy_list):
    data_iterator = ((d['K'], gsm_fidelity(d), np.sum(d['eigoccs'])) 
                     for d in iterative_data 
                     if d['L'] == L and np.isclose(d['J']/d['B'], JvB))
    K_l, fidelty_l, norms_l = zip(*sorted(data_iterator))
    fidelty_l /= np.array(norms_l)
    infidelity_l = 1 - np.array(fidelty_l)
    plt.plot(K_l, infidelity_l, 'o', label=f'iterative $J/B={JvB}$', color=c)
    
    popt, pcov = curve_fit(fitfunc, K_l[10:], infidelity_l[10:])
    plt.plot([2, 50], fitfunc(np.array([2,50]), *popt), '-', c=c, lw=1)
    plt.annotate(r'$\mathbf{\propto K^{-'+f'{popt[1]:.2f}'+r'}}$', 
                 annotate_xy, fontsize=16, zorder=10, 
                 va='center', ha='left')

    

hndlbl = [[plt.Line2D([],[],marker='',color=clrs[0]), '$J/B=.2$'],
          [plt.Line2D([],[],marker='o',ls='',color='grey'), 'iterative cooling'],
          [plt.Line2D([],[],marker='',color=clrs[1]), '$J/B=1$'],
          [plt.Line2D([],[],marker='x',ls='',color='grey'), 'reheating'],
          [plt.Line2D([],[],marker='',color=clrs[2]), '$J/B=5$'],
          [plt.Line2D([],[],marker='',ls='-',color='grey', lw=1), 'fit'],
          ]
plt.legend(*zip(*hndlbl), ncol=3, bbox_to_anchor = (0.5, 1), loc='lower center')
plt.yscale('log')
plt.xscale('log')
plt.xlabel('K')
plt.xticks([2, 10, 50], [2, 10, 50])
plt.ylabel('ground space infidelity')

#plt.savefig('../figures/Kinfidelity.pdf', bbox_inches='tight')<jupyter_output><empty_output><jupyter_text>scaling with system size L   check available data at fixed K<jupyter_code>K = 10
print(f'available data for K = {K}:')
print(' K , L, J/B ')
avail_cooling = [(d['K'], d['L'], round(d['J']/d['B'],1)) for d in cooling_data if d['K']==K]
avail_reheating = [(d['K'], d['L'], round(d['J']/d['B'],1)) for d in reheating_data if d['K']==K]
avail_iterative = [(d['K'], d['L'], round(d['J']/d['B'],1)) for d in iterative_data if d['K']==K]

from itertools import product
for K, L, JvB in np.unique(avail_cooling + avail_reheating + avail_iterative, axis=0):
    K = int(K)
    L = int(L)
    print((K, L, JvB), 
          'C' if (K, L, JvB) in avail_cooling else ' ',
          'R' if (K, L, JvB) in avail_reheating else ' ',
          'It' if (K, L, JvB) in avail_iterative else '  '
         )<jupyter_output>available data for K = 10:
 K , L, J/B 
(10, 2, 0.2) C R   
(10, 2, 1.0) C R   
(10, 2, 5.0) C R   
(10, 3, 0.2) C R   
(10, 3, 1.0) C R   
(10, 3, 5.0) C R   
(10, 4, 0.2) C R   
(10, 4, 1.0) C R   
(10, 4, 5.0) C R   
(10, 5, 0.2) C R   
(10, 5, 1.0) C R   
(10, 5, 5.0) C R   
(10, 6, 0.2) C R   
(10, 6, 1.0) C R   
(10, 6, 5.0) C R   
(10, 7, 0.2) C R It
(10, 7, 1.0) C R It
(10, 7, 5.0) C R It
(10, 8, 0.2) C R   
(10, 8, 1.0) C R   
(10, 8, 5.0) C R   
(10, 9, 1.0) C     
(10, 9, 5.0) C R<jupyter_text>energy vs L<jupyter_code>K=10
plt.figure(figsize=(9, 5))
plt.title(f'TFIM chain. Standard LogSweep(K={K}) density matrix sim')

# cooling
for JvB in [.2, 1, 5]:
    data_iterator = ((d['L'], d['energy'], np.sum(d['eigoccs'])) 
                     for d in cooling_data 
                     if d['K'] == K and np.isclose(d['J']/d['B'], JvB))
    L_l, energy_l, norms_l = zip(*sorted(data_iterator))
    energy_l /= np.array(norms_l)
    plt.plot(L_l, energy_l, 'o-', label=f'cooling $J/B={JvB}$')
    
# reheating
plt.gca().set_prop_cycle(None)
for JvB in [.2, 1, 5]:
    data_iterator = ((d['L'], d['energy'], np.sum(d['eigoccs'])) 
                     for d in reheating_data 
                     if d['K'] == K and np.isclose(d['J']/d['B'], JvB))
    L_l, energy_l, norms_l = zip(*sorted(data_iterator))
    energy_l /= np.array(norms_l)
    plt.plot(L_l, energy_l, 'x:', label=f'reheating $J/B={JvB}$')

# # iterative
# plt.gca().set_prop_cycle(None)
# for JvB in [.2, 1, 5]:
#     data_iterator = ((d['L'], d['energy'], np.sum(d['eigoccs'])) 
#                      for d in iterative_data 
#                      if d['K'] == K and np.isclose(d['J']/d['B'], JvB))
#     L_l, energy_l, norms_l = zip(*sorted(data_iterator))
#     energy_l /= np.array(norms_l)
#     plt.plot(L_l, energy_l, '+--', label=f'reheating $J/B={JvB}$')
    
plt.legend(bbox_to_anchor=(1, .5), loc='center left')
plt.ylim(top=-.95)
plt.xlabel('L')
plt.ylabel('energy')<jupyter_output><empty_output><jupyter_text>GS infidelity vs L<jupyter_code>K=10
plt.title(f'TFIM chain. Standard LogSweep(K={K}) density matrix sim')

# cooling
for JvB in [.2, 1, 5]:
    data_iterator = ((d['L'], gsm_fidelity(d), np.sum(d['eigoccs'])) 
                     for d in cooling_data 
                     if d['K'] == K and np.isclose(d['J']/d['B'], JvB))
    L_l, fidelty_l, norms_l = zip(*sorted(data_iterator))
    fidelty_l /= np.array(norms_l)
    infidelity_l = 1 - np.array(fidelty_l)
    plt.plot(L_l, infidelity_l, 'o-', label=f'cooling $J/B={JvB}$')

plt.gca().set_prop_cycle(None)

# reheating
for JvB in [.2, 1, 5]:
    data_iterator = ((d['L'], gsm_fidelity(d), np.sum(d['eigoccs'])) 
                     for d in reheating_data 
                     if d['K'] == K and np.isclose(d['J']/d['B'], JvB))
    L_l, fidelty_l, norms_l = zip(*sorted(data_iterator))
    fidelty_l /= np.array(norms_l)
    infidelity_l = 1 - np.array(fidelty_l)
    plt.plot(L_l, infidelity_l, 'x:', label=f'reheating $J/B={JvB}$')  

plt.gca().set_prop_cycle(None)

# for JvB in [.2, 1, 5]:
#     data_iterator = ((d['L'], gsm_fidelity(d), np.sum(d['eigoccs'])) 
#                      for d in iterative_data 
#                      if d['K'] == K and np.isclose(d['J']/d['B'], JvB))
#     L_l, fidelty_l, norms_l = zip(*sorted(data_iterator))
#     fidelty_l /= np.array(norms_l)
#     infidelity_l = 1 - np.array(fidelty_l)
#     plt.plot(L_l, infidelity_l, '+--', label=f'iterative $J/B={JvB}$')
    
plt.legend()
plt.yscale('log')
plt.xscale('log')
plt.xlabel('L')
plt.ylabel('ground space infidelity')<jupyter_output><empty_output><jupyter_text>test: changing K with L<jupyter_code>plt.title(f'TFIM chain. Standard LogSweep(K=L) density matrix sim')

# cooling
for JvB in [.2, 1, 5]:
    data_iterator = ((d['L'], gsm_fidelity(d), np.sum(d['eigoccs'])) 
                     for d in cooling_data 
                     if d['K'] == d['L'] and np.isclose(d['J']/d['B'], JvB))
    L_l, fidelty_l, norms_l = zip(*sorted(data_iterator))
    fidelty_l /= np.array(norms_l)
    infidelity_l = 1 - np.array(fidelty_l)
    plt.plot(L_l, infidelity_l, 'o-', label=f'cooling $J/B={JvB}$')

plt.gca().set_prop_cycle(None)

# reheating
for JvB in [.2, 1, 5]:
    data_iterator = ((d['L'], gsm_fidelity(d), np.sum(d['eigoccs'])) 
                     for d in reheating_data 
                     if d['K'] == d['L'] and np.isclose(d['J']/d['B'], JvB))
    L_l, fidelty_l, norms_l = zip(*sorted(data_iterator))
    fidelty_l /= np.array(norms_l)
    infidelity_l = 1 - np.array(fidelty_l)
    plt.plot(L_l, infidelity_l, 'x:', label=f'reheating $J/B={JvB}$')  

plt.gca().set_prop_cycle(None)

# for JvB in [.2, 1, 5]:
#     data_iterator = ((d['L'], gsm_fidelity(d), np.sum(d['eigoccs'])) 
#                      for d in iterative_data 
#                      if d['K'] == K and np.isclose(d['J']/d['B'], JvB))
#     L_l, fidelty_l, norms_l = zip(*sorted(data_iterator))
#     fidelty_l /= np.array(norms_l)
#     infidelity_l = 1 - np.array(fidelty_l)
#     plt.plot(L_l, infidelity_l, '+--', label=f'iterative $J/B={JvB}$')
    
plt.legend()
plt.yscale('log')
plt.xscale('log')
plt.xlabel('L')
plt.ylabel('ground space infidelity')<jupyter_output><empty_output><jupyter_text>Eigenstate occupation plots<jupyter_code>L = 7
JvBlist = [0.2, 1, 5]
Klist = [2, 39]

from qdclib import TFIMChain

fig, sbpl = plt.subplots(len(Klist), len(JvBlist), 
                         sharex = True, sharey = True, 
                         gridspec_kw={'hspace': 0, 'wspace': 0})

for i, JvB in enumerate(JvBlist):
    system = TFIMChain(L, JvB, 1)
    system.normalize()
    for j, K in enumerate(Klist):
        
        for d in cooling_data:
            if d['L'] == L and np.isclose(JvB, d['J']/d['B']) and d['K'] == K:
                break
        sbpl[j, i].plot(system.eigvals, d['eigoccs'], '.', label='cooling')
        
        for d in reheating_data:
            if d['L'] == L and np.isclose(JvB, d['J']/d['B']) and d['K'] == K:
                break
        sbpl[j, i].plot(system.eigvals, d['eigoccs'], '_', label='reheating')
            
        sbpl[j, i].set_yscale('log')
        sbpl[j, i].text(0.2, 0.005, f'$J/B = {JvB}$\n$K = {K}$')
        
plt.tight_layout()
sbpl[0,2].legend()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.22314293079954334, 'fraction_numerical': 0.07490884985084521, 'mean_word_length': 2.9155010814708002, 'pattern_counts': {'":': 0, '<': 88, '<?xml version=': 0, '>': 88, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '23', 'dir': 'jupyter-structured-clean-dedup', 'id': '491310', 'license': 'Apache-2.0', 'n_tokens_mistral': 14019, 'n_tokens_neox': 12485, 'n_words': 2710, 'path': 'data-analysis/TFIM-chain-logsweep-DM.ipynb', 'repo_name': 'aQaLeiden/QuantumDigitalCooling', 'size': '27153'}