Merge pull request pauxy-qmc#19 from pauxy-qmc/split_trials

Split trials

pauxy/analysis/autocorr.py

@@ -0,0 +1,70 @@
+import numpy
+import pandas as pd
+
+# Stolen from https://dfm.io/posts/autocorr/
+
+def next_pow_two(n):
+    i = 1
+    while i < n:
+        i = i << 1
+    return i
+
+def autocorr_func_1d(x, norm=True):
+    x = numpy.atleast_1d(x)
+    if len(x.shape) != 1:
+        raise ValueError("invalid dimensions for 1D autocorrelation function")
+    n = next_pow_two(len(x))
+
+    # Compute the FFT and then (from that) the auto-correlation function
+    f = numpy.fft.fft(x - numpy.mean(x), n=2*n)
+    acf = numpy.fft.ifft(f * numpy.conjugate(f))[:len(x)].real
+    acf /= 4*n
+
+    # Optionally normalize
+    if norm:
+        acf /= acf[0]
+
+    return acf
+
+# Automated windowing procedure following Sokal (1989)
+def auto_window(taus, c):
+    m = numpy.arange(len(taus)) < c * taus
+    if numpy.any(m):
+        return numpy.argmin(m)
+    return len(taus) - 1
+
+# Following the suggestion from Goodman & Weare (2010)
+def autocorr_gw2010(y, c=5.0):
+    f = autocorr_func_1d(y)
+    taus = 2.0*numpy.cumsum(f)-1.0
+    window = auto_window(taus, c)
+    return taus[window]
+
+
+def reblock_by_autocorr(y, name = "ETotal"):
+    print("# Reblock based on autocorrelation time")
+    Nmax = int(numpy.log2(len(y)))
+    Ndata = []
+    tacs = []
+    for i in range(Nmax):
+        n = int(len(y)/2**i)
+        Ndata += [n]
+        tacs += [autocorr_gw2010(y[:n])]
+
+    for n, tac in zip(Ndata, tacs):
+        print("nsamples, tac = {}, {}".format(n,tac))
+
+    block_size = int(numpy.round(numpy.max(tacs)))
+    nblocks = len(y) // block_size
+    yblocked = []
+
+    for i in range(nblocks):
+        offset = i*block_size
+        yblocked += [numpy.mean(y[offset:offset+block_size])]
+
+    yavg = numpy.mean(yblocked)
+    ystd = numpy.std(yblocked) / numpy.sqrt(nblocks)
+
+    df = pd.DataFrame({"%s_ac"%name:[yavg], "%s_error_ac"%name:[ystd], "%s_nsamp_ac"%name:[nblocks], "ac":[block_size]})
+
+    return df

pauxy/analysis/blocking.py

@@ -20,6 +20,8 @@
 from pauxy.utils.misc import get_from_dict
 from pauxy.utils.linalg import get_ortho_ao_mod
 
+from pauxy.analysis.autocorr import reblock_by_autocorr
+
 
 def average_single(frame, delete=True, multi_sym=False):
     if multi_sym:
@@ -93,7 +95,6 @@ def average_fp(frame):
                                                  2*cov_nd/(nsamp*re_num*re_den))**0.5
     return results
 
-
 def reblock_mixed(groupby, columns, verbose=False):
     analysed = []
     for group, frame in groupby:
@@ -106,6 +107,7 @@ def reblock_mixed(groupby, columns, verbose=False):
             short = short.drop(columns+drop, axis=1)
         except KeyError:
             short = short.drop(columns+['index'], axis=1)
+
         (data_len, blocked_data, covariance) = pyblock.pd_utils.reblock(short)
         reblocked = pd.DataFrame({'ETotal': [0.0]})
         for c in short.columns:
@@ -124,7 +126,14 @@ def reblock_mixed(groupby, columns, verbose=False):
             reblocked[columns[i]] = v
         analysed.append(reblocked)
 
-    return pd.concat(analysed, sort=True)
+    final = pd.concat(analysed, sort=True)
+
+    y = short["ETotal"].values
+    reblocked_ac = reblock_by_autocorr(y)
+    for c in reblocked_ac.columns:
+        final[c] = reblocked_ac[c].values
+
+    return final
 
 
 def reblock_free_projection(frame):
@@ -280,44 +289,55 @@ def analyse_back_prop(files, start_time):
         full.append(res)
     return pd.concat(full).sort_values('tau_bp')
 
-def analyse_estimates(files, start_time, multi_sim=False, verbose=False):
+def analyse_estimates(files, start_time, multi_sim=False, av_tau=False,verbose=False):
     mds = []
     basic = []
-    for f in files:
-        md = get_metadata(f)
-        read_rs = get_from_dict(md, ['psi', 'read_rs'])
-        step = get_from_dict(md, ['qmc', 'nsteps'])
-        dt = get_from_dict(md, ['qmc', 'dt'])
-        start = int(start_time/(step*dt)) + 1
-        if read_rs:
-            start = 0
-        data = extract_mixed_estimates(f, start)
-        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]
-    outfile = 'analysed_' + base
-    fmt = lambda x: "{:13.8f}".format(x)
-    print(basic_av.to_string(index=False, float_format=fmt))
-    with h5py.File(outfile, 'w') as fh5:
-        fh5['metadata'] = numpy.array(mds).astype('S')
-        try:
-            fh5['basic/estimates'] = basic_av.drop('integrals',axis=1).values.astype(float)
-        except KeyError:
-            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')
+    if av_tau:
+        data = []
+        for f in files:
+            data.append(extract_mixed_estimates(f))
+        full = pd.concat(data).groupby('Iteration')
+        av = average_tau(full)
+        print(av.apply(numpy.real).to_string())
+    else:
+        for f in files:
+            print("filename = {}".format(f))
+            md = get_metadata(f)
+            read_rs = get_from_dict(md, ['psi', 'read_rs'])
+            step = get_from_dict(md, ['qmc', 'nsteps'])
+            dt = get_from_dict(md, ['qmc', 'dt'])
+            fp = get_from_dict(md, ['propagators', 'free_projection'])
+            start = int(start_time/(step*dt)) + 1
+            if read_rs:
+                start = 0
+            data = extract_mixed_estimates(f, start)
+            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
+
+        basic = pd.concat(basic).groupby(columns)
+        if fp:
+            basic_av = reblock_free_projection(basic, columns)
+        else:
+            basic_av = reblock_mixed(basic, columns, verbose=verbose)
+
+        base = files[0].split('/')[-1]
+        outfile = 'analysed_' + base
+        fmt = lambda x: "{:13.8f}".format(x)
+        print(basic_av.to_string(index=False, float_format=fmt))
+        with h5py.File(outfile, 'w') as fh5:
+            fh5['metadata'] = numpy.array(mds).astype('S')
+            try:
+                fh5['basic/estimates'] = basic_av.drop('integrals',axis=1).values.astype(float)
+            except KeyError:
+                pass
+            fh5['basic/headers'] = numpy.array(basic_av.columns.values).astype('S')
 
 def analyse_ekt_ipea(filename, ix=None, cutoff=1e-14, screen_factor=1):
     rdm, rdm_err = average_rdm(filename, rdm_type='one_rdm', ix=ix)

pauxy/analysis/extraction.py

@@ -133,7 +133,6 @@ def get_sys_param(filename, param):
 def extract_test_data_hdf5(filename):
     """For use with testcode"""
     data = extract_mixed_estimates(filename).drop(['Iteration', 'Time'], axis=1)[::10].to_dict(orient='list')
-    # print(data)
     try:
         mrdm = extract_rdm(filename, est_type='mixed', rdm_type='one_rdm')
     except (KeyError,TypeError,AttributeError):

pauxy/estimators/ci.py

@@ -5,7 +5,7 @@
 import scipy.sparse.linalg
 import itertools
 
-def simple_fci_bose_fermi(system, nboson_max = 1, gen_dets=False, occs=None, hamil=False):
+def simple_fci_bose_fermi(system, nboson_max = 1, gen_dets=False, occs=None, hamil=False, verbose = False):
     """Very dumb FCI routine."""
     orbs = numpy.arange(system.nbasis)
 
@@ -110,42 +110,43 @@ def simple_fci_bose_fermi(system, nboson_max = 1, gen_dets=False, occs=None, ham
     Eb = eigvec[:,0].T.conj().dot(Hb.dot(eigvec[:,0]))
     Eeb = eigvec[:,0].T.conj().dot(Heb.dot(eigvec[:,0]))
 
-    for isite in range(system.nbasis):
-        rhoi = scipy.sparse.csr_matrix((ndets, ndets))
-        for i, di in enumerate(dets):
-            for d in di:
-                ii, spin_ii = map_orb(d, system.nbasis)
-                if (ii == isite):
-                    rhoi[i,i] += 1.0
-        rho = scipy.sparse.kron(Ib, rhoi)
-        nocc1 =  eigvec[:,0].T.conj().dot(rho.dot(eigvec[:,0]))
-        print("i, nocc = {}, {}".format(isite, nocc1))
-
-    for isite in range(system.nbasis):
-        bi = scipy.sparse.csr_matrix((nperms, nperms))
-        for i, iperm in enumerate(perms):
-            ni = numpy.sum(iperm)
-            offset_i = numpy.sum(blkboson[:ni+1]) # block size sum
-            if (ni == nboson_max):
-                continue
-
-            for j, jperm in enumerate(perms[offset_i:offset_i+blkboson[ni+1]]):
-                diff = numpy.array(iperm) - numpy.array(jperm)
-                ndiff = numpy.sum(numpy.abs(diff))
-                if (ndiff == 1 and diff[isite] == -1):
-                    factor = math.sqrt(numpy.array(iperm)[isite]+1)
-                    bi[i,j+offset_i] = 1.0 * factor
-
-        nib = bi.T.dot(bi)
-        ni = scipy.sparse.kron(nib, Iel)
-
-        xib = (bi + bi.T)/numpy.sqrt(2.0 * system.m * system.w0)
-        xi = scipy.sparse.kron(xib, Iel)
-
-        X =  eigvec[:,0].T.conj().dot(xi.dot(eigvec[:,0]))
-        print("i, X = {}, {}".format(isite, X))
-
-    print("# Eel, Eb, Eeb, Etot = {}, {}, {}, {}".format(Eel, Eb, Eeb, Eel+Eb+Eeb))
+    if (verbose):
+        for isite in range(system.nbasis):
+            rhoi = scipy.sparse.csr_matrix((ndets, ndets))
+            for i, di in enumerate(dets):
+                for d in di:
+                    ii, spin_ii = map_orb(d, system.nbasis)
+                    if (ii == isite):
+                        rhoi[i,i] += 1.0
+            rho = scipy.sparse.kron(Ib, rhoi)
+            nocc1 =  eigvec[:,0].T.conj().dot(rho.dot(eigvec[:,0]))
+            print("i, nocc = {}, {}".format(isite, nocc1))
+
+        for isite in range(system.nbasis):
+            bi = scipy.sparse.csr_matrix((nperms, nperms))
+            for i, iperm in enumerate(perms):
+                ni = numpy.sum(iperm)
+                offset_i = numpy.sum(blkboson[:ni+1]) # block size sum
+                if (ni == nboson_max):
+                    continue
+
+                for j, jperm in enumerate(perms[offset_i:offset_i+blkboson[ni+1]]):
+                    diff = numpy.array(iperm) - numpy.array(jperm)
+                    ndiff = numpy.sum(numpy.abs(diff))
+                    if (ndiff == 1 and diff[isite] == -1):
+                        factor = math.sqrt(numpy.array(iperm)[isite]+1)
+                        bi[i,j+offset_i] = 1.0 * factor
+
+            nib = bi.T.dot(bi)
+            ni = scipy.sparse.kron(nib, Iel)
+
+            xib = (bi + bi.T)/numpy.sqrt(2.0 * system.m * system.w0)
+            xi = scipy.sparse.kron(xib, Iel)
+
+            X =  eigvec[:,0].T.conj().dot(xi.dot(eigvec[:,0]))
+            print("i, X = {}, {}".format(isite, X))
+
+        print("# Eel, Eb, Eeb, Etot = {}, {}, {}, {}".format(Eel, Eb, Eeb, Eel+Eb+Eeb))
 
     if gen_dets:
         return (eigval, eigvec), (dets,numpy.array(oa),numpy.array(ob))

pauxy/estimators/hubbard.py

@@ -1,4 +1,53 @@
 import numpy
+def local_energy_hubbard_holstein_momentum(system, G, P, Lap, Ghalf=None):
+    r"""Calculate local energy of walker for the Hubbard-Hostein model.
+
+    Parameters
+    ----------
+    system : :class:`HubbardHolstein`
+        System information for the HubbardHolstein model.
+    G : :class:`numpy.ndarray`
+        Walker's "Green's function"
+
+    Returns
+    -------
+    (E_L(phi), T, V): tuple
+        Local, kinetic and potential energies of given walker phi.
+    """
+    # T = kinetic_lang_firsov(system.t, system.gamma_lf, P, system.nx, system.ny, system.ktwist)
+
+    Dp = numpy.array([numpy.exp(1j*system.gamma_lf*P[i]) for i in range(system.nbasis)])
+    T = numpy.zeros_like(system.T, dtype=numpy.complex128)
+    T[0] = numpy.diag(Dp).dot(system.T[0]).dot(numpy.diag(Dp.T.conj()))
+    T[1] = numpy.diag(Dp).dot(system.T[1]).dot(numpy.diag(Dp.T.conj()))
+
+    ke = numpy.sum(T[0] * G[0] + T[1] * G[1])
+
+    sqrttwomw = numpy.sqrt(2.0 * system.m * system.w0)
+    assert (system.gamma_lf * system.w0  == system.g * sqrttwomw)
+
+    Ueff = system.U + system.gamma_lf**2 * system.w0 - 2.0 * system.g * system.gamma_lf * sqrttwomw
+
+    if system.symmetric:
+        pe = -0.5*Ueff*(G[0].trace() + G[1].trace())
+
+    pe = Ueff * numpy.dot(G[0].diagonal(), G[1].diagonal())
+
+    pe_ph = - 0.5 * system.w0 ** 2 * system.m * numpy.sum(Lap)
+    ke_ph = 0.5 * numpy.sum(P*P) / system.m - 0.5 * system.w0 * system.nbasis
+
+    rho = G[0].diagonal() + G[1].diagonal()
+
+    e_eph = (system.gamma_lf**2 * system.w0 / 2.0 - system.g * system.gamma_lf * sqrttwomw) * numpy.sum(rho)
+
+    etot = ke + pe + pe_ph + ke_ph + e_eph
+
+    Eph = ke_ph + pe_ph
+    Eel = ke + pe
+    Eeb = e_eph
+
+    return (etot, ke+pe, ke_ph+pe_ph+e_eph)
+
 def local_energy_hubbard_holstein(system, G, X, Lap, Ghalf=None):
     r"""Calculate local energy of walker for the Hubbard-Hostein model.
 

pauxy/estimators/mixed.py

@@ -214,11 +214,11 @@ def update(self, system, qmc, trial, psi, step, free_projection=False):
                             E, T, V = w.local_energy(system, rchol=trial._rchol, eri=trial._eri, UVT=trial._UVT)
                         else:
                             E, T, V = 0, 0, 0
-                        self.estimates[self.names.enumer] += w.weight*E.real
+                        self.estimates[self.names.enumer] += w.weight*w.le_oratio*E.real
                         self.estimates[self.names.e1b:self.names.e2b+1] += (
-                                w.weight*numpy.array([T,V]).real
+                                w.weight*w.le_oratio*numpy.array([T,V]).real
                         )
-                        self.estimates[self.names.edenom] += w.weight
+                        self.estimates[self.names.edenom] += w.weight * w.le_oratio
                 self.estimates[self.names.uweight] += w.unscaled_weight
                 self.estimates[self.names.weight] += w.weight
                 self.estimates[self.names.ovlp] += w.weight * abs(w.ot)
@@ -442,8 +442,9 @@ def local_energy_multi_det(system, Gi, weights, two_rdm=None, rchol=None):
     denom = 0
     for w, G in zip(weights, Gi):
         # construct "local" green's functions for each component of A
-        energies += w * numpy.array(local_energy(system, G, rchol=None))
+        energies += w * numpy.array(local_energy(system, G, rchol = rchol))
         denom += w
+
     return tuple(energies/denom)
 
 def local_energy_multi_det_hh(system, Gi, weights, X, Lapi, two_rdm=None):