cleanup of library

cluster.py

@@ -9,15 +9,15 @@
 
 if __name__ == "__main__":
     cluster = False
-    job_name = 'final_h2real_'
+    job_name = 'LiH_'
     n_proc = 12
     imag_flag = False
 
     #hamiltonian = exp_loc_graph_hamiltonian(7, 1, 1)
     #hamiltonian = exp_distr_heisenberg_hamiltonian(length = 7, b_field = 5, rng_seed=1, b_rand=False)
     #hamiltonian = jellium_hamiltonian(dimensions=1, length=4, spinless=True)
-    hamiltonian = hydrogen_chain_hamiltonian(chain_length=2, bond_length=0.8)
-    #hamiltonian = LiH_hamiltonian()
+    #hamiltonian = hydrogen_chain_hamiltonian(chain_length=2, bond_length=0.8)
+    hamiltonian = LiH_hamiltonian()
     hamiltonian = normalize_hamiltonian(hamiltonian)
     norm = np.linalg.norm(np.sum(hamiltonian, axis=0), ord=2)
 
@@ -50,6 +50,7 @@
     data = {}
     data["time"] = times
     data["ham_dims"] = hamiltonian.shape
+    print(hamiltonian.shape)
     data["is_time_imaginary"] = imag_flag
     data["qdtime"] = qd_times
     data["norm"] = norm

localcluster.py

@@ -8,26 +8,30 @@
 
 if __name__ == "__main__":
     #sim params
-    cluster = True
-    n_proc = 40
-    job_name = 'local_heisenberg_9'
-    coupling = 0.0005
-    h_params = local_heisenberg_hamiltonian(length=9, coupling=coupling)
+    cluster = False
+    n_proc = 15
+    job_name = 'local_heisenberg_8_'
+    coupling = 0.00005
+    h_params = local_heisenberg_hamiltonian(length=8, coupling=coupling)
     hamiltonian = (normalize_hamiltonian(h_params[0]), h_params[1], h_params[2])
     norm = np.linalg.norm(np.sum(hamiltonian[0], axis=0), ord=2)
     print("local model shape: " + str(hamiltonian[0].shape))
 
     #time evolution
     t_i = 0.1 / norm
-    t_f= 1/norm #(3/2) * math.pi /norm
-    t_steps = 2
+    t_f= (3/2) * math.pi /norm
+    t_steps = 15
     times = list(np.geomspace(t_i, t_f, t_steps))
+    times = times[:18]
+    print(times)
     blocks = list(range(2, 5, 1))
     print(blocks)
-    block_nb = [3, 3, 3]
-    epsilon=0.00001
+    block_nb = [3, 1, 3]
+    epsilon=0.000001
+    qd_times = times[:int((3/4) * len(times))]
+    print('qd_times are : ' + str(qd_times))
 
-    partitions = ["optimal_chop", "trotter", "chop"]
+    partitions = ["trotter", "chop"]
 
     data = {}
     data['partitions'] = partitions
@@ -38,6 +42,7 @@
     data["block_nb"] = block_nb
     data["ham_dims"] = hamiltonian[0].shape
     data['coupling'] = coupling
+    data['qd_times'] = qd_times
 
     sim_params = list(product(blocks, partitions, times))
     print(sim_params)
@@ -60,38 +65,22 @@ def cost_t(params):
     for i in range(len(blocks)):
         for j in range(len(partitions)):
             data[partitions[j] + "_block_"+str(blocks[i])] = list(array_3d[i,j,:])
-
-    # #unpack results -- this will have to be modified if we include more partitions (one if statement should do it)
-    # n = len(times) # number of elements to group together
-    # count = 0
-    # index = -1
-    # b = len(blocks)
-    # for j in results:
-    #     if count % n == 0: 
-    #         index +=1 
-    #         # create a new list for every n elements
-    #         current_list = []
-    #         # associate the new list with a dictionary key
-    #         data['block ' + str(blocks[index])] = current_list
-    #     # add the current element to the current list
-    #     current_list.append(j)
-    #     count += 1
 
     def trot_cost(t):
         sim = CompositeSim(hamiltonian[0], state_rand=True, use_density_matrices=True, exact_qd=True)
         partition_sim(sim, 'trotter')
         return exact_cost(sim, t, nb=1, epsilon=epsilon)
 
-    def qd_cost(t):
-        sim = CompositeSim(hamiltonian[0], state_rand=True, use_density_matrices=True, exact_qd=True)
-        partition_sim(sim, 'qdrift')
-        return exact_cost(sim, t, nb=1, epsilon=epsilon)
+    # def qd_cost(t):
+    #     sim = CompositeSim(hamiltonian[0], state_rand=True, use_density_matrices=True, exact_qd=True)
+    #     partition_sim(sim, 'qdrift')
+    #     return exact_cost(sim, t, nb=1, epsilon=epsilon)
 
     trot_results = Parallel(n_jobs=n_proc)(delayed(trot_cost)(t) for t in times)
-    qd_results = Parallel(n_jobs=n_proc)(delayed(qd_cost)(t) for t in times[:(int(3/4 * len(times)))])
+    #qd_results = Parallel(n_jobs=n_proc)(delayed(qd_cost)(t) for t in qd_times)
 
     data['trotter'] = trot_results
-    data['qdrift'] = qd_results
+    #data['qdrift'] = qd_results
     print(data)
 
     now = datetime.now()

plotting.py

@@ -33,7 +33,7 @@ def crossover_advantage(qd_data, trot_data, comp_data, times):
     return comp_adv
 
 #plots the distribution of norms in the hamiltonian
-def ham_spec(hamiltonian_list):
+def ham_spec(hamiltonian_list, q = 75, y_line=1):
     norms = []
     index = []
     zero_norms = 0
@@ -45,10 +45,22 @@ def ham_spec(hamiltonian_list):
         if spec == 0:
             zero_norms += 1
     norms.sort()
+    q_tile = np.percentile(norms, q)
+    q_tile_loc = norms.index(q_tile)
+
+    print('q-tile is ' + str(q_tile) + ' and occurs at term '+ str(q_tile_loc))
+    print(str(len(norms)) + ' total terms in H')
+
+    trot_pot = len([i for i in norms if (i >= q_tile)])
+    print(str(trot_pot) + ' potential terms in trotter, which is ' + str(trot_pot/ len(norms) * 100) + '% of the total terms')
+
+    plt.rcParams['text.usetex'] = True
     plt.figure(figsize=(9,7))
-    plt.semilogy(index, norms, 'o-')
-    plt.xlabel("Index", size = 18)
-    plt.ylabel("Spectral Norm", size = 18)
+    plt.semilogy(index, norms, 'bs-')
+    plt.axvline(x = q_tile_loc, color = 'r', label = str(q_tile_loc) + ' Percentile')
+    plt.axhline(y = y_line, color = 'r', label = str(q_tile_loc) + ' Percentile')
+    plt.xlabel(r"$j$", size = 18)
+    plt.ylabel(r"$\left| \left| H_j \right| \right|$", size = 18)
     plt.show()
     print("There are " + str(zero_norms) + " terms with 0 spectral norm")
     return norms