Commit 6a811870 authored by Jacob August Davison's avatar Jacob August Davison
Browse files

refactored to make computations exclusively with Node objects

parent 748fe497
......@@ -10,16 +10,18 @@ from scipy.integrate import odeint, ode
import numpy as np
import time
import pickle
import tensorflow as tf
#import tensorflow as tf
#tf.enable_v2_behavior()
print("GPU available: ",tf.test.is_gpu_available())
#print("GPU available: ",tf.test.is_gpu_available())
import tracemalloc
import os, sys
#from memory_profiler import profile
import itertools
import random
sess = tf.compat.v1.Session(config=tf.compat.v1.ConfigProto(log_device_placement=True))
sess.close()
#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()
# user files
# sys.path.append('C:\\Users\\davison\\Research\\exact_diagonalization\\')
......@@ -58,12 +60,25 @@ def derivative(t, y, hamiltonian, occ_tensors, generator, flow):
E, f, G = ravel(y, hamiltonian.n_sp_states)
generator.f = f
generator.G = G
# timefi = time.time()
generator.f = tn.Node(f)
# timeff = time.time()
# timegi = time.time()
generator.G = tn.Node(G)
# timegf = time.time()
# timefli = time.time()
dE, df, dG = flow.flow(generator)
# timeflf = time.time()
# print("""
# f time: {:2.4f} s
# G time: {:2.4f} s
# Flow time: {:2.4f} s
# """.format(timeff-timefi, timegf-timegi, timeflf-timefli))
dy = unravel(dE, df, dG)
dy = unravel(dE.tensor, df.tensor, dG.tensor)
return dy
......@@ -105,6 +120,7 @@ def ravel(y, bas_len):
ravel_f = np.reshape(y[1:bas_len**2+1], (bas_len, bas_len))
ravel_G = np.reshape(y[bas_len**2+1:bas_len**2+1+bas_len**4],
(bas_len, bas_len, bas_len, bas_len))
return(ravel_E, ravel_f, ravel_G)
......@@ -119,7 +135,7 @@ def main(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0):
if ref == None:
ha = PairingHamiltonian2B(n_holes, n_particles, d=d, g=g, pb=pb)
ref = [1,1,1,1,0,0,0,0] # this is just for printing
ref = ha.reference # this is just for printing
else:
ha = PairingHamiltonian2B(n_holes, n_particles, ref=ref, d=d, g=g, pb=pb)
......@@ -146,8 +162,8 @@ def main(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0):
flowi = time.time()
# --- Solve the IM-SRG flow
y0 = unravel(ha.E, ha.f, ha.G)
y0 = unravel(ha.E.tensor, ha.f.tensor, ha.G.tensor)
#print(y0)
solver = ode(derivative,jac=None)
solver.set_integrator('vode', method='bdf', order=5, nsteps=500)
solver.set_f_params(ha, ot, wg, fl)
......@@ -216,7 +232,7 @@ if __name__ == '__main__':
# 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)
#test_exact('plots_exact_2b/', main)
main(4,4)
......
......@@ -15,7 +15,7 @@ import os, sys
#from memory_profiler import profile
import itertools
import random
import tensornetwork
import tensornetwork as tn
# user files
# sys.path.append('C:\\Users\\davison\\Research\\exact_diagonalization\\')
......@@ -258,7 +258,8 @@ def test_exact(plots_dir):
if __name__ == '__main__':
test_exact("plots3b\\")
#test_exact("plots3b/")
main3b(4,4)
#test = main3b(4,4)
#tracemalloc.start()
......
This diff is collapsed.
This diff is collapsed.
......@@ -2,7 +2,8 @@ import numpy as np
#import tensorflow as tf
# tf.enable_v2_behavior()
import tensornetwork as tn
tn.set_default_backend("tensorflow")
tn.set_default_backend("numpy")
class Hamiltonian(object):
"""Parent class for organization purposes. Ideally, all Hamiltonian
classes should inherit from this class. In this way, AssertionErrors
......@@ -238,8 +239,8 @@ class PairingHamiltonian2B(Hamiltonian):
G) -- two-body piece"""
bas1B = self.sp_basis # get the single particle basis
H1B_t = self.H1B # get the 1B tensor
H2B_t = self.H2B # get the 2B tensor
H1B_t = self.H1B.astype(np.float32) # get the 1B tensor
H2B_t = self.H2B.astype(np.float32) # get the 2B tensor
holes = self.holes # get hole states
particles = self.particles # get particle states
......@@ -248,246 +249,44 @@ class PairingHamiltonian2B(Hamiltonian):
# - Calculate 0B piece
H1B_holes = H1B_t[np.ix_(holes,holes)]
H2B_holes = H2B_t[np.ix_(holes,holes,holes,holes)]
ob_node0b = tn.Node(H1B_holes)
tb_node0b = tn.Node(H2B_holes)
ob_ii = tn.connect(ob_node0b[0],ob_node0b[1])
tb_ijij1 = tn.connect(tb_node0b[0], tb_node0b[2])
tb_ijij2 = tn.connect(tb_node0b[1], tb_node0b[3])
ob_node0b[0] ^ ob_node0b[1]
ob_contract = ob_node0b @ ob_node0b # optimized contraction
tb_node0b[0] ^ tb_node0b[2]
tb_node0b[1] ^ tb_node0b[3]
tb_contract = tb_node0b @ tb_node0b
# ob_ii = tn.connect(ob_node0b[0],ob_node0b[1])
# tb_ijij1 = tn.connect(tb_node0b[0], tb_node0b[2])
# tb_ijij2 = tn.connect(tb_node0b[1], tb_node0b[3])
flatten = tn.flatten_edges([tb_ijij1, tb_ijij2])
ob_contract = tn.contract(ob_ii).tensor.numpy()
tb_contract = 0.5*tn.contract(flatten).tensor.numpy()
# flatten = tn.flatten_edges([tb_ijij1, tb_ijij2])
# ob_contract = tn.contract(ob_ii).tensor#.tensor.numpy()
# tb_contract = 0.5*tn.contract(flatten).tensor#.tensor.numpy()
E = ob_contract + tb_contract
E = E.astype(np.float32)
E = ob_contract + tn.Node(0.5)*tb_contract
#E = E.astype(np.float32)
# - Calculate 1B piece
ob_node1b = tn.Node(H1B_t)
tb_node1b = tn.Node(H2B_t[np.ix_(bas1B,holes,bas1B,holes)])
tb_ihjh = tn.connect(tb_node1b[1], tb_node1b[3])
tb_contract = tn.contract(tb_ihjh)
tb_node1b[1] ^ tb_node1b[3]
tb_contract = tb_node1b @ tb_node1b
f = ob_node1b.tensor.numpy() + tb_contract.tensor.numpy()
f = f.astype(np.float32)
G = H2B_t
# tb_ihjh = tn.connect(tb_node1b[1], tb_node1b[3])
# tb_contract = tn.contract(tb_ihjh)
#f = ob_node1b.tensor.numpy() + tb_contract.tensor.numpy()
f = (ob_node1b + tb_contract)
#f = f.astype(np.float32)
# - Calculate 2B piece
G = tn.Node(H2B_t)
return (E, f, G)
# class PairingHamiltonian3B(PairingHamiltonian2B):
#
# def __init__(self, n_hole_states, n_particle_states, d=1.0, g=0.5, pb=0.0):
# """Class constructor. Instantiate PairingHamiltonian3B object.
#
# Arguments:
#
# n_hole_states -- number of holes states in the single particle basis
# n_particle_states -- number of particles states in the single particle basis
#
# Keyword arguments:
#
# d -- the energy level spacing (default: 1.0)
# g -- the pairing strength (default: 0.5)
# pb -- strength of the pair-breaking term (operates in double particle basis) (default: 0.0)"""
#
# self._d = d
# self._g = g
# self._pb = pb
# self._reference = np.append(np.ones(n_hole_states), np.zeros(n_particle_states))
# self._holes = np.arange(n_hole_states, dtype=np.int64)
#
# self._n_sp_states = n_hole_states + n_particle_states
# self._particles = np.arange(n_hole_states,self.n_sp_states, dtype=np.int64)
# self._sp_basis = np.append(self.holes, self.particles)
#
# self._H1B, self._H2B, self._H3B = self.construct()
# self._E, self._f, self._G = self.normal_order()
#
# @property
# def H3B(self):
# """Returns:
#
# H2B -- three-body (rank 6) tensor defined by construct()."""
# return self._H3B
#
# def delta2B(self, p,q,r,s):
# """Inherited from parent."""
# return super().delta2B(p,q,r,s)
#
# def deltaPB(self, p,q,r,s):
# """Inherited from parent."""
# return super().deltaPB(p,q,r,s)
#
# def delta3B(self, p,q,r,s,t,u):
# """Determines if a two-body tensor elements should be zero,
# positive, or negative. This behavior is dictated by the three-
# body excitation term in pairing Hamiltonian.
#
# Arguments:
#
# p,q,r,s,t,u -- indices in single-particle basis"""
#
# pp = np.floor_divide(p,2)
# qp = np.floor_divide(q,2)
# rp = np.floor_divide(r,2)
# sp = np.floor_divide(s,2)
# tp = np.floor_divide(t,2)
# up = np.floor_divide(u,2)
#
# ps = 1 if p%2==0 else -1
# qs = 1 if q%2==0 else -1
# rs = 1 if r%2==0 else -1
# ss = 1 if s%2==0 else -1
# ts = 1 if t%2==0 else -1
# us = 1 if u%2==0 else -1
#
# lhs_p = np.array([pp, qp, rp])
# lhs_s = np.array([ps, qs, rs])
# rhs_p = np.array([sp, tp, up])
# rhs_s = np.array([ss, ts, us])
#
# p_equiv = np.equal(lhs_p, rhs_p)
# s_equiv = np.equal(lhs_s, rhs_s)
# p_where = np.where(p_equiv==True)[0]
# s_where = np.where(s_equiv==True)[0]
# # print(p_where)
# result = 0
# result_str ='0'
# if len(p_where) == 2 and (lhs_p[p_where[0]] == lhs_p[p_where[1]] and rhs_p[p_where[0]] == rhs_p[p_where[1]]):
# # print("in")
# result = 1
#
# p_where_n = np.where(p_equiv==False)[0]
# l_equiv = np.equal(lhs_p[p_where_n], lhs_p[p_where])
# r_equiv = np.equal(rhs_p[p_where_n], rhs_p[p_where])
# # print(l_equiv, r_equiv)
# if(l_equiv.any() or r_equiv.any()):
# result *= 0
#
# l_spins = lhs_s[p_where]
# r_spins = rhs_s[p_where]
# # print(l_spins, r_spins)
#
# if l_spins[0] == l_spins[1] or r_spins[0] == r_spins[1]:
# result *= 0
# result_str += '0'
# if l_spins[0] == r_spins[0] and l_spins[1] == r_spins[1]:
# result *= 1
# result_str += '1'
# if l_spins[0] == r_spins[1] and l_spins[1] == r_spins[0]:
# result *= -1
# result_str += '-1'
#
# if lhs_p[p_where_n] != rhs_p[p_where_n] and lhs_s[p_where_n] != rhs_p[p_where_n]:
# result *= 1
# result_str += '1'
# else:
# result *= 0
# result_str += '0'
#
# return result
#
# def construct(self):
# """Constructs the one- and two-body pieces of the pairing
# Hamiltonian.
#
# Returns:
#
# (H1B, -- one-body tensor elements (defined by one-body operator)
# H2B) -- two-body tensor elements (defined by two-body operator)"""
#
# bas1B = self.sp_basis # get the single particle basis
#
# # one body part of Hamiltonian is floor-division of basis index
# # matrix elements are (P-1) where P is energy level
# H1B = np.diag(self.d*np.floor_divide(bas1B,2))
#
# # two body part of Hamiltonian constructed from four indices
# # with non-zero elements defined by pairing term
# H2B = np.zeros(np.ones(4, dtype=np.int64)*self.n_sp_states)
# for p in bas1B:
# for q in bas1B:
# for r in bas1B:
# for s in bas1B:
# H2B[p,q,r,s] += self.g*0.5*self.delta2B(p,q,r,s)
# H2B[p,q,r,s] += self.pb*0.5*self.deltaPB(p,q,r,s)
#
# W=1
# H3B = np.zeros(np.ones(6, dtype=np.int64)*self.n_sp_states)
# for p in bas1B:
# for q in bas1B:
# for r in bas1B:
# for s in bas1B:
# for t in bas1B:
# for u in bas1B:
# H3B[p,q,r,s,t,u] += W*(1./6.)*self.delta3B(p,q,r,s,t,u)
# # print(H3B[0,1,2,3,4,5])
# return (H1B, H2B, H3B)
# # if lhs_p[p_where_n] != rhs_p[p_where_n]:
# # if np.array_equal(s_where, [])
#
# def normal_order(self):
# """Normal-orders the pairing Hamiltonian.
#
# Returns:
#
# (E, -- zero-body piece
# f, -- one-body piece
# G) -- two-body piece"""
#
# bas1B = self.sp_basis # get the single particle basis
# H1B_t = self.H1B # get the 1B tensor
# H2B_t = self.H2B # get the 2B tensor
# H3B_t = self.H3B
# holes = self.holes # get hole states
# particles = self.particles # get particle states
#
# net = TensorNetwork()
#
# # - Calculate 0B piece
# H1B_holes = H1B_t[np.ix_(holes,holes)]
# H2B_holes = H2B_t[np.ix_(holes,holes,holes,holes)]
# H3B_holes = H3B_t[np.ix_(holes,holes,holes,holes,holes,holes)]
#
# ob_node0b = net.add_node(H1B_holes)
# tb_node0b = net.add_node(H2B_holes)
#
# ob_ii = net.connect(ob_node0b[0],ob_node0b[1])
# tb_ijij1 = net.connect(tb_node0b[0], tb_node0b[2])
# tb_ijij2 = net.connect(tb_node0b[1], tb_node0b[3])
#
# flatten = net.flatten_edges([tb_ijij1, tb_ijij2])
# ob_contract = net.contract(ob_ii).tensor.numpy()
# tb_contract = 0.5*net.contract(flatten).tensor.numpy()
#
# E = ob_contract + tb_contract
#
#
# # - Calculate 1B piece
# ob_node1b = net.add_node(H1B_t)
# tb_node1b = net.add_node(H2B_t[np.ix_(bas1B,holes,bas1B,holes)])
#
# tb_ihjh = net.connect(tb_node1b[1], tb_node1b[3])
# tb_contract = net.contract(tb_ihjh)
#
# f = ob_node1b.tensor.numpy() + tb_contract.tensor.numpy()
#
# G = H2B_t
#
# del net
#
# return (E, f, G)
#
# # test = PairingHamiltonian3B(4,4)
# # bas1B = test.sp_basis
# # for p in bas1B:
# # for q in bas1B:
# # for r in bas1B:
# # for s in bas1B:
# # for t in bas1B:
# # for u in bas1B:
# # if test.H3B[p,q,r,s,t,u] != 0:
# # print('{:0.4f} | {:d} {:d} {:d} {:d} {:d} {:d}'.format(test.H3B[p,q,r,s,t,u], p,q,r,s,t,u))
return (E, f, G)
import numpy as np
import os
from oop_imsrg.hamiltonian import *
import sys
import numba
class OccupationTensors(object):
"""Functions as a container for important occupation tensors
......@@ -31,6 +35,11 @@ class OccupationTensors(object):
self._occI = self.__get_occI()
self._occJ = self.__get_occJ()
if not os.path.exists("occ_storage/"):
os.mkdir("occ_storage/")
@property
def occA(self):
"""Returns:
......@@ -112,8 +121,9 @@ class OccupationTensors(object):
return self._occJ
# ---- BUILD OCCUPATION TENSORS ---
# ---- BUILD OCCUPATION TENSORS ---
#@jit#(nopython=True)
def __get_occA(self, flag=0):
"""Builds the occupation tensor occA.
......@@ -130,11 +140,48 @@ class OccupationTensors(object):
n = len(bas1B)
if flag == 0: # default
occA = np.zeros((n,n,n,n),dtype=np.float32)
for a in bas1B:
for b in bas1B:
occA[a,b,a,b] = ref[a] - ref[b]
# occA = np.zeros((n,n,n,n),dtype=np.float32)
# for a in bas1B:
# for b in bas1B:
# occA[a,b,a,b] = ref[a] - ref[b]
# print(sys.getsizeof(occA)/10**6)
# TENSOR TRAIN DECOMPOSITION
# We find a TT-decomposition by hand, because the the tensor we want
# to decompose is small enough to do so. This is an exact decomposition
# of the rank 2 tensor described by n_a-n_b.
#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)))
Gab = tn.ncon([Ga,Gb], [(-1,1),(1,-2)])
final = tn.outer_product(Gab, tn.Node(np.ones((n,n)))).tensor
# PARALLELIZE NESTED LOOPS FOR BETTER PERFORMANCE
@numba.jit(nopython=True)
def enforce_delta(n, tensor):
bas1B = range(n)
for a in bas1B:
for b in bas1B:
for c in bas1B:
for d in bas1B:
if not(a == c and b == d):
tensor[a,b,c,d] = 0
return tensor
occA = tn.Node(enforce_delta(n, final))
# print(sys.getsizeof(occA)/10**6)
if flag == 1:
occA = np.zeros((n,n),dtype=np.float32)
......@@ -161,11 +208,41 @@ class OccupationTensors(object):
n = len(bas1B)
if flag == 0: # default
occB = np.zeros((n,n,n,n),dtype=np.float32)
# occB = np.zeros((n,n,n,n),dtype=np.float32)
for a in bas1B:
for b in bas1B:
occB[a,b,a,b] = 1 - ref[a] - ref[b]
# for a in bas1B:
# for b in bas1B:
# occB[a,b,a,b] = 1 - ref[a] - ref[b]
#occB = tn.Node(1 - self.occA.tensor)
# TENSOR TRAIN DECOMPOSITION of rank 4 tensor with elements given
# by 1 - n_a - n_b.
#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)))
Gab = tn.ncon([Ga,Gb], [(-1,1),(1,-2)])
final = tn.outer_product(Gab, tn.Node(np.ones((n,n)))).tensor
# PARALLELIZE NESTED LOOPS FOR BETTER PERFORMANCE
@numba.jit(nopython=True)
def enforce_delta(n, tensor):
bas1B = range(n)
for a in bas1B:
for b in bas1B:
for c in bas1B:
for d in bas1B:
if not(a == c and b == d):
tensor[a,b,c,d] = 0
return tensor
occB = tn.Node(enforce_delta(n, final))
if flag == 1:
occB = np.zeros((n,n),dtype=np.float32)
......@@ -191,27 +268,66 @@ 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)
for a in bas1B:
for b in bas1B:
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]
# occC = np.zeros((n,n,n,n,n,n),dtype=np.float32)
# for a in bas1B:
# for b in bas1B:
# 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]
# print(sys.getsizeof(occC)/10**6)
# TENSOR TRAIN DECOMPOSITION of rank 6 tensor with elements
# given by n_a*n_b + (1 - n_b - n_a)*n_c.
#Ga = tn.Node(np.array([ [1,1,1,1], [1,1,1,1], [1,1,1,1], [1,1,1,1], [0,1,1,0], [0,1,1,0], [0,1,1,0], [0,1,1,0] ]))
# Gb = tn.Node(np.array([ [[1,0,0,0], [0,1,0,0],[0,0,-1,0],[0,0,0,-1]],
# [[1,0,0,0], [0,1,0,0],[0,0,-1,0],[0,0,0,-1]],
# [[1,0,0,0], [0,1,0,0],[0,0,-1,0],[0,0,0,-1]],
# [[1,0,0,0], [0,1,0,0],[0,0,-1,0],[0,0,0,-1]],
# [[0,0,0,0], [0,1,0,0],[0,0,0,0],[0,0,0,-1]],
# [[0,0,0,0], [0,1,0,0],[0,0,0,0],[0,0,0,-1]],
# [[0,0,0,0], [0,1,0,0],[0,0,0,0],[0,0,0,-1]],
# [[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)
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)
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)
Gc = tn.Node(np.transpose(Gc1))
Gabc = tn.ncon([Ga, Gb, Gc], [(-1, 1), (-2, 1, 2), (2, -3)])
final = tn.outer_product(Gabc, tn.Node(np.ones((n,n,n)))).tensor
@numba.jit(nopython=True)#, parallel=True)
def enforce_delta(n, tensor):
bas1B = range(n)
for a in bas1B:
for b in bas1B:
for c in bas1B:
for d in bas1B:
for e in bas1B:
for f in bas1B:
if not(a == d and b == e and c == f):
tensor[a,b,c,d,e,f] = 0
return tensor