Skip to content

Commit

Permalink
Merge pull request pauxy-qmc#15 from pauxy-qmc/sri
Browse files Browse the repository at this point in the history 
SRI migration from pauxy-dev
  • Loading branch information
@fdmalone
fdmalone authored Aug 10, 2020
2 parents 0be3d12 + 52018d6 commit 559c6f4
Show file tree
Hide file tree
Showing 16 changed files with 422 additions and 216 deletions.
13 changes: 11 additions & 2 deletions pauxy/analysis/blocking.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -124,7 +124,6 @@ def reblock_mixed(groupby, columns, verbose=False):
reblocked[columns[i]] = v
analysed.append(reblocked)


return pd.concat(analysed, sort=True)


Expand Down Expand Up @@ -296,6 +295,15 @@ def analyse_estimates(files, start_time, multi_sim=False, verbose=False):
columns = set_info(data, md)
basic.append(data.drop('Iteration', axis=1))
mds.append(md)
new_columns = []
for c in columns:
if (c != "E_T"):
new_columns += [c]
columns = new_columns
if (len(files) > 1):
print("multi simulations detected")
print("grouping based on everything other than E_T")
print("columns = {}".format(columns))
basic = pd.concat(basic).groupby(columns)
basic_av = reblock_mixed(basic, columns, verbose=verbose)
base = files[0].split('/')[-1]
Expand All @@ -307,7 +315,8 @@ def analyse_estimates(files, start_time, multi_sim=False, verbose=False):
try:
fh5['basic/estimates'] = basic_av.drop('integrals',axis=1).values.astype(float)
except KeyError:
pass
print("integrals does not exist under the problem class")
fh5['basic/estimates'] = basic_av.values.astype(float)
fh5['basic/headers'] = numpy.array(basic_av.columns.values).astype('S')

def analyse_ekt_ipea(filename, ix=None, cutoff=1e-14, screen_factor=1):
Expand Down
37 changes: 22 additions & 15 deletions pauxy/analysis/extraction.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,19 +13,23 @@ def extract_data_sets(files, group, estimator, raw=False):

def extract_data(filename, group, estimator, raw=False):
fp = get_param(filename, ['propagators', 'free_projection'])
with h5py.File(filename, 'r') as fh5:
dsets = list(fh5[group][estimator].keys())
data = numpy.array([fh5[group][estimator][d][:] for d in dsets])
if 'rdm' in estimator or raw:
return data
else:
header = fh5[group]['headers'][:]
header = numpy.array([h.decode('utf-8') for h in header])
df = pd.DataFrame(data)
df.columns = header
if not fp:
df = df.apply(numpy.real)
return df
try:
with h5py.File(filename, 'r') as fh5:
dsets = list(fh5[group][estimator].keys())
data = numpy.array([fh5[group][estimator][d][:] for d in dsets])
if 'rdm' in estimator or raw:
return data
else:
header = fh5[group]['headers'][:]
header = numpy.array([h.decode('utf-8') for h in header])
df = pd.DataFrame(data)
df.columns = header
if not fp:
df = df.apply(numpy.real)
return df
except:
print("Problem with {}".format(filename))
exit()

def extract_mixed_estimates(filename, skip=0):
return extract_data(filename, 'basic', 'energies')[skip:]
Expand Down Expand Up @@ -110,8 +114,11 @@ def set_info(frame, md):
return list(frame.columns[ncols:])

def get_metadata(filename):
with h5py.File(filename, 'r') as fh5:
metadata = json.loads(fh5['metadata'][()])
try:
with h5py.File(filename, 'r') as fh5:
metadata = json.loads(fh5['metadata'][()])
except:
print("# problem with file = {}".format(filename))
return metadata

def get_param(filename, param):
Expand Down
225 changes: 218 additions & 7 deletions pauxy/estimators/generic.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,17 +32,25 @@ def local_energy_generic(h1e, eri, G, ecore=0.0, Ghalf=None):
return (e1+e2+ecore, e1+ecore, e2)

def local_energy_generic_opt(system, G, Ghalf=None):
# import cProfile
# pr = cProfile.Profile()
# pr.enable()

# Element wise multiplication.
e1b = numpy.sum(system.H1[0]*G[0]) + numpy.sum(system.H1[1]*G[1])
Gup = Ghalf[0].ravel()
Gdn = Ghalf[1].ravel()
euu = 0.5 * Gup.dot(system.vakbl[0].dot(Gup))
edd = 0.5 * Gdn.dot(system.vakbl[1].dot(Gdn))
# TODO: Fix this. Very dumb.
eos = numpy.dot((system.rchol_vecs[0].T).dot(Gup),
(system.rchol_vecs[1].T).dot(Gdn))

# eos = numpy.dot((system.rchol_vecs[0].T).dot(Gup),
# (system.rchol_vecs[1].T).dot(Gdn))
eos = Gup.dot(system.vakbl[2].dot(Gdn))
e2b = euu + edd + eos #eud + edu

# pr.disable()
# pr.print_stats(sort='tottime')

return (e1b + e2b + system.ecore, e1b + system.ecore, e2b)

def local_energy_generic_cholesky_opt(system, G, Ghalf, rchol):
Expand All @@ -62,10 +70,17 @@ def local_energy_generic_cholesky_opt(system, G, Ghalf, rchol):
(E, T, V): tuple
Local, kinetic and potential energies.
"""
#import cProfile
#pr = cProfile.Profile()
#pr.enable()

# Element wise multiplication.
e1b = numpy.sum(system.H1[0]*G[0]) + numpy.sum(system.H1[1]*G[1])
nalpha, nbeta = system.nup, system.ndown
nbasis = system.nbasis
if rchol is not None:
naux = rchol.shape[1]

Ga, Gb = Ghalf[0], Ghalf[1]
Xa = rchol[:nalpha*nbasis].T.dot(Ga.ravel())
Xb = rchol[nalpha*nbasis:].T.dot(Gb.ravel())
Expand All @@ -75,12 +90,208 @@ def local_energy_generic_cholesky_opt(system, G, Ghalf, rchol):
rchol_a, rchol_b = rchol[:nalpha*nbasis], rchol[nalpha*nbasis:]
# T_{abn} = \sum_k Theta_{ak} LL_{ak,n}
# LL_{ak,n} = \sum_i L_{ik,n} A^*_{ia}
Ta = numpy.tensordot(Ga, rchol_a.reshape((nalpha,nbasis,-1)), axes=((1),(1)))
exxa = numpy.tensordot(Ta, Ta, axes=((0,1,2),(1,0,2)))
Tb = numpy.tensordot(Gb, rchol_b.reshape((nbeta,nbasis,-1)), axes=((1),(1)))
exxb = numpy.tensordot(Tb, Tb, axes=((0,1,2),(1,0,2)))
# rchol_a = rchol_a.reshape((nalpha,nbasis, naux))
# rchol_b = rchol_b.reshape((nbeta,nbasis, naux))

# Ta = numpy.tensordot(Ga, rchol_a, axes=((1),(1)))
# exxa = numpy.tensordot(Ta, Ta, axes=((0,1,2),(1,0,2)))
# Tb = numpy.tensordot(Gb, rchol_b, axes=((1),(1)))
# exxb = numpy.tensordot(Tb, Tb, axes=((0,1,2),(1,0,2)))

rchol_a = rchol_a.T
rchol_b = rchol_b.T
Ta = numpy.zeros((naux, nalpha, nalpha), dtype=rchol_a.dtype)
Tb = numpy.zeros((naux, nbeta, nbeta), dtype=rchol_b.dtype)
GaT = Ga.T
GbT = Gb.T
for x in range(naux):
rmi_a = rchol_a[x].reshape((nalpha,nbasis))
Ta[x] = rmi_a.dot(GaT)
rmi_b = rchol_b[x].reshape((nbeta,nbasis))
Tb[x] = rmi_b.dot(GbT)
exxa = numpy.tensordot(Ta, Ta, axes=((0,1,2),(0,2,1)))
exxb = numpy.tensordot(Tb, Tb, axes=((0,1,2),(0,2,1)))

exx = exxa + exxb
e2b = 0.5 * (ecoul - exx)

#pr.disable()
#pr.print_stats(sort='tottime')
return (e1b + e2b + system.ecore, e1b + system.ecore, e2b)

# def local_energy_generic_cholesky_opt_stochastic(system, G, nsamples, Ghalf=None, rchol=None):
# r"""Calculate local for generic two-body hamiltonian.

# This uses the cholesky decomposed two-electron integrals.

# Parameters
# ----------
# system : :class:`hubbard`
# System information for the hubbard model.
# G : :class:`numpy.ndarray`
# Walker's "green's function"

# Returns
# -------
# (E, T, V): tuple
# Local, kinetic and potential energies.
# """

# # Element wise multiplication.
# e1b = numpy.sum(system.H1[0]*G[0]) + numpy.sum(system.H1[1]*G[1])

# if rchol is None:
# rchol = system.rchol_vecs

# nalpha, nbeta= system.nup, system.ndown
# nbasis = system.nbasis

# Ga, Gb = Ghalf[0], Ghalf[1]
# Xa = rchol[0].T.dot(Ga.ravel())
# Xb = rchol[1].T.dot(Gb.ravel())
# ecoul = numpy.dot(Xa,Xa)
# ecoul += numpy.dot(Xb,Xb)
# ecoul += 2*numpy.dot(Xa,Xb)

# # O(ON s)
# # Xa = numpy.tensordot(ra, Ga, axes=((0,1),(0,1)))
# # Xb = numpy.tensordot(rb, Gb, axes=((0,1),(0,1)))
# # ecoul = numpy.dot(Xa,Xa)
# # ecoul += numpy.dot(Xb,Xb)
# # ecoul += 2*numpy.dot(Xa,Xb)

# naux = rchol[0].shape[-1]
# theta = numpy.zeros((naux,nsamples), dtype=numpy.int64)
# for i in range(nsamples):
# theta[:,i] = (2*numpy.random.randint(0,2,size=(naux))-1)

# if system.sparse:
# rchol_a, rchol_b = [rchol[0].toarray(), rchol[1].toarray()]
# else:
# rchol_a, rchol_b = rchol[0], rchol[1]

# rchol_a = rchol_a.reshape((nalpha,nbasis, naux))
# rchol_b = rchol_b.reshape((nbeta,nbasis, naux))

# # ra = numpy.tensordot(rchol_a, theta, axes=((2),(0))) * numpy.sqrt(1.0/nsamples)
# # rb = numpy.tensordot(rchol_b, theta, axes=((2),(0))) * numpy.sqrt(1.0/nsamples)
# ra = numpy.einsum("ipX,Xs->ips",rchol_a, theta, optimize=True) * numpy.sqrt(1.0/nsamples)
# rb = numpy.einsum("ipX,Xs->ips",rchol_b, theta, optimize=True) * numpy.sqrt(1.0/nsamples)
# Gra = numpy.einsum("kq,lqx->lkx", Ga, ra, optimize=True)
# Grb = numpy.einsum("kq,lqx->lkx", Gb, rb, optimize=True)

# # Gra = numpy.tensordot(Ga, ra, axes=((1),(1))).transpose(1,0,2)
# # Grb = numpy.tensordot(Gb, rb, axes=((1),(1))).transpose(1,0,2)
# exxa = numpy.tensordot(Gra, Gra, axes=((0,1,2),(1,0,2)))
# exxb = numpy.tensordot(Grb, Grb, axes=((0,1,2),(1,0,2)))

# exx = exxa + exxb
# e2b = 0.5 * (ecoul - exx)

# return (e1b + e2b + system.ecore, e1b + system.ecore, e2b)
def local_energy_generic_cholesky_opt_stochastic(system, G, nsamples, Ghalf=None, rchol=None, C0=None,\
ecoul0 = None, exxa0 = None, exxb0 = None):
r"""Calculate local for generic two-body hamiltonian.
This uses the cholesky decomposed two-electron integrals.
Parameters
----------
system : :class:`hubbard`
System information for the hubbard model.
G : :class:`numpy.ndarray`
Walker's "green's function"
Returns
-------
(E, T, V): tuple
Local, kinetic and potential energies.
"""
#import cProfile
#pr = cProfile.Profile()
#pr.enable()

if (type(C0) == numpy.ndarray):
control = True
else:
control = False

# Element wise multiplication.
e1b = numpy.sum(system.H1[0]*G[0]) + numpy.sum(system.H1[1]*G[1])
if rchol is None:
rchol = system.rchol_vecs
nalpha, nbeta= system.nup, system.ndown
nbasis = system.nbasis
Ga, Gb = Ghalf[0], Ghalf[1]
Xa = rchol[0].T.dot(Ga.ravel())
Xb = rchol[1].T.dot(Gb.ravel())
ecoul = numpy.dot(Xa,Xa)
ecoul += numpy.dot(Xb,Xb)
ecoul += 2*numpy.dot(Xa,Xb)
if system.sparse:
rchol_a, rchol_b = [rchol[0].toarray(), rchol[1].toarray()]
else:
rchol_a, rchol_b = rchol[0], rchol[1]

# T_{abn} = \sum_k Theta_{ak} LL_{ak,n}
# LL_{ak,n} = \sum_i L_{ik,n} A^*_{ia}

naux = rchol_a.shape[-1]

theta = numpy.zeros((naux,nsamples), dtype=numpy.int64)
for i in range(nsamples):
theta[:,i] = (2*numpy.random.randint(0,2,size=(naux))-1)

if (control):

ra = rchol_a.dot(theta).T * numpy.sqrt(1.0/nsamples)
rb = rchol_b.dot(theta).T * numpy.sqrt(1.0/nsamples)

Ta0 = numpy.zeros((nsamples, nalpha, nalpha), dtype=rchol_a.dtype)
Tb0 = numpy.zeros((nsamples, nbeta, nbeta), dtype=rchol_b.dtype)

Ta = numpy.zeros((nsamples, nalpha, nalpha), dtype=rchol_a.dtype)
Tb = numpy.zeros((nsamples, nbeta, nbeta), dtype=rchol_b.dtype)

G0aT = C0[:,:system.nup]
G0bT = C0[:,system.nup:]

GaT = Ga.T
GbT = Gb.T

for x in range(nsamples):
rmi_a = ra[x].reshape((nalpha,nbasis))
rmi_b = rb[x].reshape((nbeta,nbasis))

Ta0[x] = rmi_a.dot(G0aT)
Tb0[x] = rmi_b.dot(G0bT)
Ta[x] = rmi_a.dot(GaT)
Tb[x] = rmi_b.dot(GbT)

exxa_hf = numpy.tensordot(Ta0, Ta0, axes=((0,1,2),(0,2,1)))
exxb_hf = numpy.tensordot(Tb0, Tb0, axes=((0,1,2),(0,2,1)))

exxa_corr = numpy.tensordot(Ta, Ta, axes=((0,1,2),(0,2,1)))
exxb_corr = numpy.tensordot(Tb, Tb, axes=((0,1,2),(0,2,1)))

exxa = exxa0 + (exxa_corr - exxa_hf)
exxb = exxb0 + (exxb_corr - exxb_hf)

else:
rchol_a = rchol_a.reshape((nalpha,nbasis, naux))
rchol_b = rchol_b.reshape((nbeta,nbasis, naux))

ra = numpy.einsum("ipX,Xs->ips",rchol_a, theta, optimize=True) * numpy.sqrt(1.0/nsamples)
Gra = numpy.einsum("kq,lqx->lkx", Ga, ra, optimize=True)
exxa = numpy.tensordot(Gra, Gra, axes=((0,1,2),(1,0,2)))

rb = numpy.einsum("ipX,Xs->ips",rchol_b, theta, optimize=True) * numpy.sqrt(1.0/nsamples)
Grb = numpy.einsum("kq,lqx->lkx", Gb, rb, optimize=True)
exxb = numpy.tensordot(Grb, Grb, axes=((0,1,2),(1,0,2)))

exx = exxa + exxb
e2b = 0.5 * (ecoul - exx)

#pr.disable()
#pr.print_stats(sort='tottime')

return (e1b + e2b + system.ecore, e1b + system.ecore, e2b)

def local_energy_generic_cholesky(system, G, Ghalf=None):
Expand Down
Loading

0 comments on commit 559c6f4

Please sign in to comment.