memory optimization and removed for loop in occE

main.py

@@ -20,6 +20,7 @@ import itertools
 import random
 sess = tf.compat.v1.Session(config=tf.compat.v1.ConfigProto(log_device_placement=True))
 sess.close()
+
 # user files
 # sys.path.append('C:\\Users\\davison\\Research\\exact_diagonalization\\')
 from oop_imsrg.hamiltonian import *
@@ -195,38 +196,53 @@ def main(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0):
     
     return (convergence, iters, d, g, pb, n_holes+n_particles, s_vals, E_vals, time_str)
 
-def exact_diagonalization(d, g):
-    """Result of exact diagonalization in spin=0 block of
-    pairing Hamiltonian, given 8 single particle states (4 hole states
-    and 4 particle states).
-
-    Arguments:
-
-    d -- energy level spacing
-    g -- pairing strength
-
-    Returns:
-
-    E -- ground state energy
-    """
-    H = [[2*d-g, -g/2, -g/2, -g/2, -g/2, 0],
-         [-g/2, 4*d-g, -g/2, -g/2, 0, -g/2],
-         [-g/2, -g/2, 6*d-g, 0, -g/2, -g/2],
-         [-g/2, -g/2, 0, 6*d-g, -g/2, -g/2],
-         [-g/2, 0, -g/2, -g/2, 8*d-g, -g/2],
-         [0, -g/2, -g/2, -g/2, -g/2, 10*d-g]]
-
-    w, v = np.linalg.eig(H)
-    E = w[0]
-
-    return E
-
 
 if __name__ == '__main__':
 
     #test_refs('logs_refs\\')
-    print(ci_matrix.exact_diagonalization(1.0, 0.5,0.0))
-    test = main(4,4)
-
-#     # test_exact('plots_exact\\')
-    # print(exact_diagonalization(1.0,0.5))
+    # print(ci_matrix.exact_diagonalization(1.0, 0.5,0.0))
+   #test = main(4,4)
+    # occt = OccupationTensors(np.array([0,1,2,3,4,5,6,7]), np.array([1,1,1,1,0,0,0,0]))
+    # bas1B = [0,1,2,3,4,5,6,7]
+    # occE = occt.occE
+    # for a in bas1B:
+    #     for b in bas1B:
+    #         for c in bas1B:
+    #             for d in bas1B:
+    #                 for e in bas1B:
+    #                     for f in bas1B:
+    #                         val = occE[a,b,c,d,e,f]
+    #                         if val != 0:
+    #                             with open('occE_nonzero.txt', 'a') as fi:
+    #                                 fi.write('%s %s %s %s %s %s -- %s\n' % (a,b,c,d,e,f,val))
+
+    #test_exact('plots_exact_2b/', main)
+    
+    bas1B = np.array([0,1,2,3,4,5,6,7])
+    ref = np.array([1,1,1,1,0,0,0,0],dtype=np.float16)
+    occt = OccupationTensors(bas1B,ref)
+    occF = occt.occF
+    
+    test1 = np.einsum('i,j,k,l,m->ijklm',ref,ref,(1-ref),(1-ref),(1-ref))
+    test1e = np.einsum('ijklm,nopqr->ijklmnopqr',test1,test1)
+    test2 = np.einsum('i,j,k,l,m->ijklm',(1-ref),(1-ref),ref,ref,ref)
+    test2e = np.einsum('ijklm,nopqr->ijklmnopqr',test2,test2)
+    test = test1e + test2e
+    # test = np.einsum('ijklm,nopqr->ijklmnopqr',test3,test3)
+    print(occF.dtype, test.dtype)
+    print(occF.shape, test.shape)
+    print(np.array_equal(occF, test))
+    
+    # for a in bas1B:
+    #     for b in bas1B:
+    #         for c in bas1B:
+    #             for d in bas1B:
+    #                 for e in bas1B:
+    #                     for f in bas1B:
+    #                         val1 = test[a,b,c,d,e,f]
+    #                         val2 = occC[a,b,c,d,e,f]
+    #                         if val1 != val2:
+    #                             print('conflict found for ',a,b,c,d,e,f'-',val1,val2)
+                            
+    #                         with open('occE_nonzero_test.txt', 'a') as fi:
+    #                             fi.write('%s %s %s %s %s -- %s\n' % (a,b,c,d,e,val))

main3b.py

@@ -156,7 +156,7 @@ def main3b(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0):
 
     # --- Solve the IM-SRG flow
     y0 = unravel(ha.E, ha.f, ha.G, wg.W)
-    print('sum(Ws**2)',np.sum(wg.W**2))
+    
     solver = ode(derivative,jac=None)
     solver.set_integrator('vode', method='bdf', order=5, nsteps=500)
     solver.set_f_params(ha, ot, wg, fl)
@@ -175,9 +175,10 @@ def main3b(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0):
         Es, fs, Gs, Ws = ravel(ys, ha.n_sp_states)
         s_vals.append(solver.t)
         E_vals.append(Es)
-        print('sum(Ws**2)',np.sum(Ws**2))
+        
         iters += 1
-
+        # if iters == 1:
+        #     break
         if iters %10 == 0: print("iter: {:>6d} \t scale param: {:0.4f} \t E = {:0.9f}".format(iters, solver.t, Es))
 
         if len(E_vals) > 100 and abs(E_vals[-1] - E_vals[-2]) < 10**-8 and E_vals[-1] != E_vals[0]:
@@ -219,47 +220,46 @@ def test_exact(plots_dir):
 
     g_vals = np.linspace(start, stop, num)
 
-    for pb in g_vals:
-        E_corrs = []
-        E_exacts = []
-        for g in g_vals:
-            data = main3b(4,4, d=1.0, g=g, pb=0.0)
-            E_vals = data[7]
-            E_corr = E_vals[-1]
-            E_exact = exact_diagonalization(1.0, g)
-
-            E_corrs.append(E_corr - (2-g))
-            E_exacts.append(E_exact - (2-g))
-
-            plt.figure(figsize=[12,8])
-            plt.plot(data[6], data[7])
-            plt.ylabel('Energy')
-            plt.xlabel('scale parameter')
-            plt.title('Convergence for \n g={:2.4f}, pb={:2.4f}'.format(g,pb))
+    E_corrs = []
+    E_exacts = []
+    for g in g_vals:
+        data = main3b(4,4, d=1.0, g=g, pb=0.0)
+        E_vals = data[7]
+        E_corr = E_vals[-1]
+        E_exact = ci_matrix.exact_diagonalization(1.0, g, 0.0)
 
-            pb_plots_dir = plots_dir+'pb{:2.4f}\\'.format(pb)
-            if not os.path.exists(pb_plots_dir):
-                os.mkdir(pb_plots_dir)
-
-            plt.savefig(pb_plots_dir+'g{:2.4f}_pb{:2.4f}.png'.format(g,pb))
-            plt.close()
+        E_corrs.append(E_corr - (2-g))
+        E_exacts.append(E_exact - (2-g))
 
         plt.figure(figsize=[12,8])
-        plt.plot(g_vals, E_exacts, marker='s')
-        plt.plot(g_vals, E_corrs, marker='v')
-        plt.ylabel('E$_{corr}$')
-        plt.xlabel('g')
-        plt.legend(['exact', 'IMSRG(3)'])
-        plt.title('Correlation energy with pb = {:2.4f}'.format(pb))
-        plt.savefig(plots_dir+'pb{:2.4f}.png'.format(pb))
+        plt.plot(data[6], data[7])
+        plt.ylabel('Energy')
+        plt.xlabel('scale parameter')
+        plt.title('Convergence for \n g={:2.4f}, pb={:2.4f}'.format(g,0.0))
+        
+        pb_plots_dir = plots_dir+'pb{:2.4f}\\'.format(0.0)
+        if not os.path.exists(pb_plots_dir):
+            os.mkdir(pb_plots_dir)
+
+        plt.savefig(pb_plots_dir+'g{:2.4f}_pb{:2.4f}.png'.format(g,0.0))
         plt.close()
-        print(E_exacts)
-        break
+
+    plt.figure(figsize=[12,8])
+    plt.plot(g_vals, E_exacts, marker='s')
+    plt.plot(g_vals, E_corrs, marker='v')
+    plt.ylabel('E$_{corr}$')
+    plt.xlabel('g')
+    plt.legend(['exact', 'IMSRG(3)'])
+    plt.title('Correlation energy with pb = {:2.4f}'.format(pb))
+    plt.savefig(plots_dir+'pb{:2.4f}.png'.format(pb))
+    plt.close()
+    print(E_exacts)
+
 
 
 if __name__ == '__main__':
-    # test_exact("plots3b\\")
-    test = main3b(4,4)
+    test_exact("plots3b\\")
+    #test = main3b(4,4)
     
     #tracemalloc.start()
     

oop_imsrg/hamiltonian.py

@@ -36,9 +36,10 @@ class PairingHamiltonian2B(Hamiltonian):
         self._g = g
         self._pb = pb
         if ref == None:
-            self._reference = np.append(np.ones(n_hole_states), np.zeros(n_particle_states))
+            self._reference = np.append(np.ones(n_hole_states,dtype=np.float32),
+                                        np.zeros(n_particle_states,dtype=np.float32))
         else:
-            self._reference = ref
+            self._reference = np.asarray(ref,dtype=np.float32)
 
         self._holes = np.arange(n_hole_states, dtype=np.int32)
         self._n_sp_states = n_hole_states + n_particle_states

oop_imsrg/occupation_tensors.py

@@ -191,6 +191,18 @@ class OccupationTensors(object):
         ref = self._reference
         n = len(bas1B)
 
+        # occC1 = np.einsum('i,j,k->ijk', ref,ref,(1-ref))
+        # occC2 = np.einsum('i,j,k->ijk', (1-ref),(1-ref),ref)
+        # occC3 = occC1 + occC2
+        
+        # if flag == 0: # default
+        #     occC = np.einsum('ijk,lmn->ijklmn', occC3, occC3)
+        #     return occC
+        
+        # if flag == 1:
+        #     occC = occC3
+        #     return occC
+
         if flag == 0: # default
             occC = np.zeros((n,n,n,n,n,n),dtype=np.float32)
 
@@ -199,8 +211,7 @@ class OccupationTensors(object):
                     for c in bas1B:
                         occC[a,b,c,a,b,c] = ref[a]*ref[b]*(1-ref[c]) + \
                                             (1-ref[a])*(1-ref[b])*ref[c]
-
-
+        
         if flag == 1:
             occC = np.zeros((n,n,n),dtype=np.float32)
 
@@ -225,7 +236,14 @@ class OccupationTensors(object):
         bas1B = self._sp_basis
         ref = self._reference
         n = len(bas1B)
-
+        # occD = np.einsum('i,j,k,l->ijkl',ref,ref,(1-ref),(1-ref))
+        
+        # if flag == 0: # default
+        #     occD = np.einsum('ijkl,mnop->ijklmnop', occD, occD)
+        #     return occD
+        # if flag == 1:
+        #     return occD
+        
         if flag == 0: # default
             occD = np.zeros((n,n,n,n,n,n,n,n),dtype=np.float32)
 
@@ -245,9 +263,8 @@ class OccupationTensors(object):
                         for d in bas1B:
                             occD[a,b,c,d] = ref[a]*ref[b]*\
                                             (1-ref[c])*(1-ref[d])
-
         return occD
-
+                            
 # ---- ALL ABOVE REQUIRED FOR IMSRG(2) ---
 
     def __get_occE(self):
@@ -260,19 +277,22 @@ class OccupationTensors(object):
         ref = self._reference
         n = len(bas1B)
 
-        occE = np.zeros((n,n,n,n,n,n),dtype=np.float32)
+        occE = np.einsum('i,j,k,l,m,n->ijklmn',ref,ref,ref,(1-ref),(1-ref),(1-ref))
+        
+        # occE = np.zeros((n,n,n,n,n,n),dtype=np.float32)
 
-        for a in bas1B:
-            for b in bas1B:
-                for c in bas1B:
-                    for d in bas1B:
-                        for e in bas1B:
-                            for f in bas1B:
-                                occE[a,b,c,d,e,f] = ref[a]*ref[b]*ref[c]*\
-                                                    (1-ref[d])*(1-ref[e])*\
-                                                    (1-ref[f])
+        # for a in bas1B:
+        #     for b in bas1B:
+        #         for c in bas1B:
+        #             for d in bas1B:
+        #                 for e in bas1B:
+        #                     for f in bas1B:
+        #                         occE[a,b,c,d,e,f] = ref[a]*ref[b]*ref[c]*\
+        #                                             (1-ref[d])*(1-ref[e])*\
+        #                                             (1-ref[f])
 
         return occE
+    
     def __get_occF(self):
         """Builds the occupation tensor __get_occF. Treat as a rank 10 tensor.
 
@@ -284,6 +304,11 @@ class OccupationTensors(object):
         ref = self._reference
         n = len(bas1B)
 
+        # occF1 = np.einsum('i,j,k,l,m->ijklm',ref,ref,(1-ref),(1-ref),(1-ref))
+        # occF2 = np.einsum('i,j,k,l,m->ijklm',(1-ref),(1-ref),ref,ref,ref)
+        # occF3 = occF1 + occF2
+        # occF = np.einsum('ijklm,nopqr->ijklmnopqr',occF3,occF3)
+        
         occF = np.zeros((n,n,n,n,n,n,n,n,n,n),dtype=np.float32)
 
         for a in bas1B: