changes to occ tensors; ensemble ref support

benchmarking_wd/imsrg_pairing.py

@@ -23,6 +23,7 @@ from scipy.integrate import odeint, ode
 from sys import argv
 
 import time
+import pickle
 #-----------------------------------------------------------------------------------
 # basis and index functions
 #-----------------------------------------------------------------------------------
@@ -898,7 +899,7 @@ def main(n_holes, g=0.5):
 
   # set up initial Hamiltonian
   H1B, H2B = pairing_hamiltonian(delta, g, user_data)
-  
+  pickle.dump((H1B, H2B), open('comparison.p','wb'))
   E, f, Gamma = normal_order(H1B, H2B, user_data) 
   #print(E)
 
@@ -953,12 +954,12 @@ def main(n_holes, g=0.5):
 # make executable
 #------------------------------------------------------------------------------
 if __name__ == "__main__": 
-  for n in range(2,14,2):
-    ti = time.time()
-    main(n)
-    tf = time.time()
-    print('Total time: {:2.5f}'.format(tf-ti))
-
+  # for n in range(2,14,2):
+  #   ti = time.time()
+  #   main(n)
+  #   tf = time.time()
+  #   print('Total time: {:2.5f}'.format(tf-ti))
+  main(4)
 import pytest
 
 def test02ph_run(benchmark):

main.py

@@ -10,6 +10,8 @@ from scipy.integrate import odeint, ode
 import numpy as np
 import time
 import pickle
+
+import matplotlib.pyplot as plt
 #import tensorflow as tf
 #tf.enable_v2_behavior()
 #print("GPU available: ",tf.test.is_gpu_available())
@@ -18,7 +20,7 @@ import os, sys
 #from memory_profiler import profile
 import itertools
 import random
-#import tensornetwork as tn
+import tensornetwork as tn
 #tn.set_default_backend("tensorflow")
 #sess = tf.compat.v1.Session(config=tf.compat.v1.ConfigProto(log_device_placement=True))
 #sess.close()
@@ -34,6 +36,20 @@ from oop_imsrg.plot_data import *
 import oop_imsrg.ci_pairing.cipy_pairing_plus_ph as ci_matrix
 from oop_imsrg.tests2B import *
 
+def get_vacuum_coeffs(E, f, G, basis, holes):
+
+    H2B = G
+    H1B = f - np.trace(G[np.ix_(basis,holes,basis,holes)], axis1=1,axis2=3) 
+
+    Gnode = tn.Node(G[np.ix_(holes,holes,holes,holes)])
+    Gnode[0] ^ Gnode[2]
+    Gnode[1] ^ Gnode[3]
+    result = Gnode @ Gnode
+    H0B = E - np.trace(f[np.ix_(holes,holes)]) + 0.5*result.tensor
+
+    return (H0B, H1B, H2B)
+
+
 # @profile
 def derivative(t, y, hamiltonian, occ_tensors, generator, flow):
     """Defines the derivative to pass into ode object.
@@ -125,7 +141,7 @@ def ravel(y, bas_len):
     return(ravel_E, ravel_f, ravel_G)
 
 # @profile
-def main(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0, verbose=1):
+def main(n_holes, n_particles, ref=[], d=1.0, g=0.5, pb=0.0, verbose=1):
     """Main method uses scipy.integrate.ode to solve the IMSRG(2) flow
     equations."""
 
@@ -133,7 +149,7 @@ def main(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0, verbose=1):
 
     initi = time.time() # start instantiation timer
 
-    if ref == None:
+    if ref == []:
         ha = PairingHamiltonian2B(n_holes, n_particles, d=d, g=g, pb=pb)
         ref = ha.reference # this is just for printing
     else:
@@ -165,6 +181,8 @@ def main(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0, verbose=1):
 
     # --- Solve the IM-SRG flow
     y0 = unravel(ha.E, ha.f, ha.G)
+    coeffs = get_vacuum_coeffs(ha.E, ha.f, ha.G, ha.sp_basis, ha.holes)
+    pickle.dump( coeffs, open( "./vac_coeffs_unevolved.p", "wb" ) )
 
     solver = ode(derivative,jac=None)
     solver.set_integrator('vode', method='bdf', order=5, nsteps=500)
@@ -190,10 +208,17 @@ def main(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0, verbose=1):
         #     break
         if iters %10 == 0 and verbose: print("iter: {:>6d} \t scale param: {:0.4f} \t E = {:0.8f}".format(iters, solver.t, Es))
 
+#        if iters %20 == 0 and verbose:
+#            coeffs = get_vacuum_coeffs(Es, fs, Gs, ha.sp_basis, ha.holes)
+#            pickle.dump( coeffs, open( "mixed_state_test/pickled_coeffs/vac_coeffs_s{}.p".format(iters), "wb" ) )
+
         if len(E_vals) > 100 and abs(E_vals[-1] - E_vals[-2]) < 10**-8 and E_vals[-1] != E_vals[0]:
 
             if verbose: print("---- Energy converged at iter {:>06d} with energy {:1.8f}\n".format(iters,E_vals[-1]))
             convergence = 1
+            coeffs = get_vacuum_coeffs(Es, fs, Gs, ha.sp_basis, ha.holes)
+            #pickle.dump( coeffs, open( "mixed_state_test/pickled_coeffs/vac_coeffs_evolved.p", "wb" ) )
+            pickle.dump(coeffs, open('vac_coeffs_evolved.p', 'wb'))
             break
 
         if len(E_vals) > 100 and abs(E_vals[-1] - E_vals[-2]) > 1:
@@ -206,6 +231,8 @@ def main(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0, verbose=1):
 
         if iters % 1000 == 0 and verbose:
             print('Iteration {:>06d}'.format(iters))
+
+        
     flowf = time.time()
     end = time.time()
     time_str = "{:2.5f}".format(end-start)
@@ -221,9 +248,91 @@ if __name__ == '__main__':
 
 
     #test_exact('plots_exact_2b/', main)
+    # refs = list(map("".join, itertools.permutations('11110000')))
+    # refs = list(dict.fromkeys(refs)) # remove duplicates
+    # refs = [list(map(int, list(ref))) for ref in refs]
+
+
+    # TESTING MIXED STATE --------------------------------------------
+    # refs = [[1,1,1,1,0,0,0,0],[1,1,0,0,1,1,0,0],[1,1,0,0,0,0,1,1],
+    #         [0,0,1,1,1,1,0,0],[0,0,1,1,0,0,1,1],[0,0,0,0,1,1,1,1]]
+    
+
+    # gsws = [0.99, 0.985,0.975,0.95,0.9,0.85,0.8]#,0.75,0.7,0.65,0.6,0.55,0.5]
+    
+    # E_data_full = []
+    # for gsw in gsws:
+    #     E_data = []
+    #     count = 0
+    #     rsw = (1-gsw)
+    #     for i in range(len(refs)):
+    #         ref = np.zeros_like(refs[0])
+    #         E = 0.0
+    #         if i == 0: 
+    #             ref = refs[0]
+    #             data = main(4,4,ref=ref,verbose=0)
+    #             E = (data[7])[-1]
+    #             count += 1
+    #         else:
+    #             esum = np.zeros_like(refs[0])
+    #             for state in refs[1:count+1]:
+    #                 esum = np.add(esum,state)
+    #             print(esum)
+    #             ref = gsw*np.array(refs[0]) + rsw/count*(esum)
+    #             data = main(4,4,ref=ref,verbose=0)
+    #             E = (data[7])[-1]
+    #             count += 1
+    #         print("{:0.8f}, {}, {f}".format(E, sum(ref), f=ref))
+    #         E_data.append(E)
+    #     E_data_full.append(E_data)
+
+    # exact = ci_matrix.exact_diagonalization(1.0, 0.5, 0.0)
+
+    # pickle.dump( E_data_full, open( "save.p", "wb" ) )
+
+    # ----------------------------------------------------------------
+
+    refs = [[1,1,1,1,0,0,0,0],[1,1,0,0,1,1,0,0],[1,1,0,0,0,0,1,1],
+            [0,0,1,1,1,1,0,0],[0,0,1,1,0,0,1,1]]
+
+    ref = 0.985*np.asarray(refs[0]) + (1.0-0.985)/4.0*(np.asarray(refs[1]) + np.asarray(refs[2]) + np.asarray(refs[3]) + np.asarray(refs[4]))
+    main(4,4, g=0.95, ref=[1,1,1,1,0,0,0,0])
+    # H1B_true, H2B_true = pickle.load(open('comparison.p','rb'))
+    # H1B, H2B = pickle.load(open('vac_coeffs_unevolved.p', 'rb'))
+    # print(H1B, H1B_true)
+
+    
+    #print(np.array_equal(H2B_true, H2B))
+    #main(4,4)
+    
+
+    # ref1 = 0.9*refs[0] + 0.1*refs[1]
+    # data = main(4,4,ref=ref1)
+    # E1 = data[7]
+
+    # ref2 = 0.9*refs[0] + 0.05*(refs[1] + refs[2])
+    # data = main(4,4,ref=ref2)
+        
+#    ref = 0.9*np.array([1,1,1,1,0,0,0,0])+0.1*np.array([0,1,1,1,1,0,0,0])
+    
+    #main(4,4,ref=ref)
+
+    # exact = ci_matrix.exact_diagonalization(1.0, 0.5, 0.0)
+
+    # fig = plt.figure()
+
+    # for i in range(0,10):
+    #     data = main(4,4)
+    #     s_vals = data[6]
+    #     E_vals = data[7]/exact
+    #     plt.plot(s_vals, E_vals)
+    # plt.ylim([1+10**-8, 1+10**-6])
+    # plt.xlabel('scale param (s)')
+    # plt.ylabel('energy')
+    # plt.show()
 
-    main(4,4)
 
+    #-- TESTING TIMINGS -------------------------
     #holes = int(sys.argv[1])
 
     # print('convergence,iters,d,g,pb,num states,GS energy,total time')

oop_imsrg/generator.py

 #import tensorflow as tf
 # tf.enable_v2_behavior()
 #import tensornetwork as tn
-#import numpy as np
+import numpy as np
 from oop_imsrg.hamiltonian import *
 from oop_imsrg.occupation_tensors import *
 #tn.set_default_backend("tensorflow") 
@@ -42,7 +42,13 @@ class WegnerGenerator(Generator):
         self._occC = occ_t.occC
         self._occD = occ_t.occD
 
-#        self._occRef1 = 
+        ref = h.reference
+        n = len(self._holes)+len(self._particles)
+        Ga = tn.Node(np.transpose(np.append(ref[np.newaxis,:], np.zeros((1,n)), axis=0).astype(float)))
+        Gb = tn.Node(np.append(ref[::-1][np.newaxis,:],np.zeros((1,n)), axis=0).astype(float))
+        
+        self._occRef1 = tn.ncon([Ga,Gb], [(-1,1),(1,-2)])                                                                    # n_a(1-n_b)
+        self._occRef2 = tn.ncon([tn.Node(np.transpose(Gb.tensor)), tn.Node(np.transpose(Ga.tensor))], [(-1,1),(1,-2)])       # (1-n_a)n_b
 
     @property
     def f(self):
@@ -86,16 +92,26 @@ class WegnerGenerator(Generator):
         holes = self._holes
         particles = self._particles
 
+        occRef1 = self._occRef1
+        occRef2 = self._occRef2
+
+        occD = self._occD
+        occDt = tn.Node(np.transpose(occD.tensor,[2,3,0,1]))
+
         # - Decouple off-diagonal 1B and 2B pieces
  #       fod1 = tn.ncon([])
         fod = np.zeros(f.shape, dtype=np.float32)
-        fod[np.ix_(particles, holes)] += f[np.ix_(particles, holes)]
-        fod[np.ix_(holes, particles)] += f[np.ix_(holes, particles)]
+        fod += np.multiply(occRef2.tensor, f)
+        fod += np.multiply(occRef1.tensor, f)
+        # fod[np.ix_(particles, holes)] += f[np.ix_(particles, holes)]
+        # fod[np.ix_(holes, particles)] += f[np.ix_(holes, particles)]
         fd = f - fod
 
         God = np.zeros(G.shape, dtype=np.float32)
-        God[np.ix_(particles, particles, holes, holes)] += G[np.ix_(particles, particles, holes, holes)]
-        God[np.ix_(holes, holes, particles, particles)] += G[np.ix_(holes, holes, particles, particles)]
+        God += np.multiply(occDt.tensor, G)
+        God += np.multiply(occD.tensor, G)
+        # God[np.ix_(particles, particles, holes, holes)] += G[np.ix_(particles, particles, holes, holes)]
+        # God[np.ix_(holes, holes, particles, particles)] += G[np.ix_(holes, holes, particles, particles)]
         Gd = G - God
 
         return (fd, fod, Gd, God)
@@ -143,7 +159,7 @@ class WegnerGenerator(Generator):
         sum2_1b_1 = tn.ncon([fdPrime, God], [(1,2), (2,-1,1,-2)])
         sum2_1b_2 = tn.ncon([fodPrime, Gd], [(1,2), (2,-1,1,-2)])
         sum2_1b = sum2_1b_1 - sum2_1b_2
-
+        
         # sum2_1b_1 = tn.ncon([fd, God], [(0, 1), (1, -1, 0, -2)])#.numpy()
         # sum2_1b_2 = tn.ncon([fod, Gd], [(0, 1), (1, -1, 0, -2)])#.numpy()
         # sum2_1b_3 = sum2_1b_1 - sum2_1b_2
@@ -152,7 +168,7 @@ class WegnerGenerator(Generator):
         # third term
         #sum3_1b_1 = tn.ncon([Gd, occC, God], [(6,-1,4,5), (4,5,6,1,2,3), (1,2,3,-2)])#.numpy()
         #sum3_1b = sum3_1b_1 - np.transpose(sum3_1b_1)
-
+        
         sum3_1b_1 = np.multiply(occC.tensor, God)#np.multiply(tn.outer_product(tn.Node(occC), tn.Node(np.ones(8))).tensor, God)
         sum3_1b_2 = tn.ncon([Gd, God], [(3,-1,1,2),(1,2,3,-2)])
         sum3_1b = sum3_1b_2 - np.transpose(sum3_1b_2)

oop_imsrg/hamiltonian.py

@@ -18,7 +18,7 @@ class Hamiltonian(object):
 class PairingHamiltonian2B(Hamiltonian):
     """Generate the two-body pairing Hamiltonian. Inherits from Hamiltonian."""
 
-    def __init__(self, n_hole_states, n_particle_states, ref=None, d=1.0, g=0.5, pb=0.0):
+    def __init__(self, n_hole_states, n_particle_states, ref=[], d=1.0, g=0.5, pb=0.0):
         """Class constructor. Instantiate PairingHamiltonian2B object.
 
         Arguments:
@@ -37,7 +37,7 @@ class PairingHamiltonian2B(Hamiltonian):
         self._d = d
         self._g = g
         self._pb = pb
-        if ref == None:
+        if ref == []:
             self._reference = np.append(np.ones(n_hole_states,dtype=np.float32),
                                         np.zeros(n_particle_states,dtype=np.float32))
         else:

oop_imsrg/occupation_tensors.py

@@ -12,6 +12,8 @@ class OccupationTensors(object):
     index is occupied, or 0 if the same index is unoccupied,
     in the reference state."""
 
+    DATA_TYPE = float
+
     def __init__(self, sp_basis, reference):
         """Class constructor. Instantiates an OccupationTensors object.
 
@@ -184,8 +186,8 @@ class OccupationTensors(object):
             #Ga = tn.Node(np.array([ [1,1], [1,1], [1,1], [1,1], [0,1], [0,1], [0,1], [0,1] ]))
             #Gb = tn.Node(np.transpose(np.array([ [1,-1], [1,-1], [1,-1], [1,-1], [1,0], [1,0], [1,0], [1,0] ])))
 
-            Ga = tn.Node(np.append(ref[:,np.newaxis], np.ones((n,1)),axis=1).astype(int))
-            Gb = tn.Node(np.transpose(np.append(np.ones((n,1)), -1*ref[:,np.newaxis],axis=1).astype(int)))
+            Ga = tn.Node(np.append(ref[:,np.newaxis], np.ones((n,1)),axis=1).astype(self.DATA_TYPE))
+            Gb = tn.Node(np.transpose(np.append(np.ones((n,1)), -1*ref[:,np.newaxis],axis=1).astype(self.DATA_TYPE)))
             Gab = tn.ncon([Ga,Gb], [(-1,1),(1,-2)])
             #final = tn.outer_product(Gab, tn.Node(np.ones((n,n)))).tensor
 
@@ -216,8 +218,8 @@ class OccupationTensors(object):
             #     for b in bas1B:
             #         occA[a,b] = ref[a] - ref[b]
 
-            Ga = tn.Node(np.append(ref[:,np.newaxis], np.ones((n,1)),axis=1).astype(int))
-            Gb = tn.Node(np.transpose(np.append(np.ones((n,1)), -1*ref[:,np.newaxis],axis=1).astype(int)))
+            Ga = tn.Node(np.append(ref[:,np.newaxis], np.ones((n,1)),axis=1).astype(self.DATA_TYPE))
+            Gb = tn.Node(np.transpose(np.append(np.ones((n,1)), -1*ref[:,np.newaxis],axis=1).astype(self.DATA_TYPE)))
             Gab = tn.ncon([Ga,Gb], [(-1,1),(1,-2)])
             
             final = tn.outer_product(Gab, tn.Node(np.ones((n,n))))
@@ -256,8 +258,8 @@ class OccupationTensors(object):
             #Ga = tn.Node(np.array([ [0,1], [0,1], [0,1], [0,1], [1,1], [1,1], [1,1], [1,1] ]))
             #Gb = tn.Node(np.transpose(np.array([ [1,-1], [1,-1], [1,-1], [1,-1], [1,0], [1,0], [1,0], [1,0] ])))
 
-            Ga = tn.Node(np.append(1-ref[:,np.newaxis], np.ones((n,1)),axis=1).astype(int))
-            Gb = tn.Node(np.transpose(np.append(np.ones((n,1)), -1*ref[:,np.newaxis],axis=1).astype(int)))
+            Ga = tn.Node(np.append(1-ref[:,np.newaxis], np.ones((n,1)),axis=1).astype(self.DATA_TYPE))
+            Gb = tn.Node(np.transpose(np.append(np.ones((n,1)), -1*ref[:,np.newaxis],axis=1).astype(self.DATA_TYPE)))
             Gab = tn.ncon([Ga,Gb], [(-1,1),(1,-2)])
             #final = tn.outer_product(Gab, tn.Node(np.ones((n,n)))).tensor
 
@@ -286,8 +288,8 @@ class OccupationTensors(object):
             #     for b in bas1B:
             #         occB[a,b] = 1 - ref[a] - ref[b]
 
-            Ga = tn.Node(np.append(1-ref[:,np.newaxis], np.ones((n,1)),axis=1).astype(int))
-            Gb = tn.Node(np.transpose(np.append(np.ones((n,1)), -1*ref[:,np.newaxis],axis=1).astype(int)))
+            Ga = tn.Node(np.append(1-ref[:,np.newaxis], np.ones((n,1)),axis=1).astype(self.DATA_TYPE))
+            Gb = tn.Node(np.transpose(np.append(np.ones((n,1)), -1*ref[:,np.newaxis],axis=1).astype(self.DATA_TYPE)))
             Gab = tn.ncon([Ga,Gb], [(-1,1),(1,-2)])
             final = tn.outer_product(Gab, tn.Node(np.ones((n,n))))
 
@@ -336,16 +338,16 @@ class OccupationTensors(object):
             #                         [[0,0,0,0], [0,1,0,0],[0,0,0,0],[0,0,0,-1]] ]))
             #Gc = tn.Node(np.transpose(np.array([ [1,1,1,1], [1,1,1,1], [1,1,1,1], [1,1,1,1], [1,0,0,0], [1,0,0,0], [1,0,0,0], [1,0,0,0] ])))
 
-            Ga1 = np.append(ref[:,np.newaxis], np.ones((n,2)),axis=1).astype(int)
-            Ga2 = np.append(Ga1, ref[:,np.newaxis],axis=1).astype(int)
+            Ga1 = np.append(ref[:,np.newaxis], np.ones((n,2)),axis=1).astype(self.DATA_TYPE)
+            Ga2 = np.append(Ga1, ref[:,np.newaxis],axis=1).astype(self.DATA_TYPE)
             Ga = tn.Node(Ga2) 
 
-            Gb1= np.append(ref[np.newaxis,np.newaxis,:],np.zeros((1,1,n)),axis=1).astype(int)
-            Gb2= np.append(Gb1, np.append(np.zeros((1,1,n)), np.ones((1,1,n)),axis=1), axis=0).astype(int)
+            Gb1= np.append(ref[np.newaxis,np.newaxis,:],np.zeros((1,1,n)),axis=1).astype(self.DATA_TYPE)
+            Gb2= np.append(Gb1, np.append(np.zeros((1,1,n)), np.ones((1,1,n)),axis=1), axis=0).astype(self.DATA_TYPE)
             Gb3 = np.array([[1,0],[0,-1]])
             Gb = tn.Node(np.kron(Gb3, np.transpose(Gb2)))
            
-            Gc1 = np.append(np.ones((n,1)), np.repeat(ref[:,np.newaxis],3,axis=1), axis=1).astype(int)
+            Gc1 = np.append(np.ones((n,1)), np.repeat(ref[:,np.newaxis],3,axis=1), axis=1).astype(self.DATA_TYPE)
             Gc = tn.Node(np.transpose(Gc1))
 
             Gabc = tn.ncon([Ga, Gb, Gc], [(-1, 1), (-2, 1, 2), (2, -3)])
@@ -419,9 +421,9 @@ class OccupationTensors(object):
             # Gc = tn.Node(np.array([ [0,0],[0,0],[0,0],[0,0],[1,0],[1,0],[1,0],[1,0] ]))
             # Gd = tn.Node(np.transpose(np.array([ [0,0],[0,0],[0,0],[0,0],[1,0],[1,0],[1,0],[1,0] ])))
 
-            Ga = tn.Node(np.transpose(np.append(ref[np.newaxis,:], np.zeros((1,n)), axis=0).astype(int)))
+            Ga = tn.Node(np.transpose(np.append(ref[np.newaxis,:], np.zeros((1,n)), axis=0).astype(self.DATA_TYPE)))
             Gb = tn.Node(np.transpose(Ga.tensor))
-            Gc = tn.Node(np.transpose(np.append(ref[::-1][np.newaxis,:],np.zeros((1,n)), axis=0).astype(int)))
+            Gc = tn.Node(np.transpose(np.append(ref[::-1][np.newaxis,:],np.zeros((1,n)), axis=0).astype(self.DATA_TYPE)))
             Gd = tn.Node(np.transpose(Gc.tensor))
 
             Gabcd = tn.ncon([Ga,Gb,Gc,Gd], [(-1,1),(1,-2),(-3,2),(2,-4)])
@@ -463,9 +465,9 @@ class OccupationTensors(object):
             # Gc = tn.Node(np.array([ [0,0],[0,0],[0,0],[0,0],[1,0],[1,0],[1,0],[1,0] ]))
             # Gd = tn.Node(np.transpose(np.array([ [0,0],[0,0],[0,0],[0,0],[1,0],[1,0],[1,0],[1,0] ])))
 
-            Ga = tn.Node(np.transpose(np.append(ref[np.newaxis,:], np.zeros((1,n)), axis=0).astype(int)))
+            Ga = tn.Node(np.transpose(np.append(ref[np.newaxis,:], np.zeros((1,n)), axis=0).astype(self.DATA_TYPE)))
             Gb = tn.Node(np.transpose(Ga.tensor))
-            Gc = tn.Node(np.transpose(np.append(ref[::-1][np.newaxis,:],np.zeros((1,n)), axis=0).astype(int)))
+            Gc = tn.Node(np.transpose(np.append(ref[::-1][np.newaxis,:],np.zeros((1,n)), axis=0).astype(self.DATA_TYPE)))
             Gd = tn.Node(np.transpose(Gc.tensor))
 
             Gabcd = tn.ncon([Ga,Gb,Gc,Gd], [(-1,1),(1,-2),(-3,2),(2,-4)])