mirror of https://github.com/QMCPACK/qmcpack.git
97 lines
2.4 KiB
Python
Executable File
97 lines
2.4 KiB
Python
Executable File
#!/usr/bin/env python
|
|
from numpy import *
|
|
import sys
|
|
|
|
def corr(i,x,mean,var):
|
|
N=len(x)
|
|
if var==0:#if the variance is 0 return an effectively infinity corr
|
|
return 1e100
|
|
corr=1.0/var*1.0/(N-i)*sum((x[0:N-i]-mean)*(x[i:N]-mean))
|
|
return corr
|
|
|
|
def Stats(x):
|
|
N=len(x)
|
|
mean=sum(x)/(N+0.0)
|
|
xSquared=x*x
|
|
var=sum(xSquared)/(N+0.0)-mean*mean
|
|
i=0
|
|
tempC=0.5
|
|
kappa=0.0
|
|
while (tempC>0 and i<(N-1)):
|
|
kappa=kappa+2.0*tempC
|
|
i=i+1
|
|
tempC=corr(i,x,mean,var)
|
|
if kappa == 0.0:
|
|
kappa = 1.0
|
|
Neff=(N+0.0)/(kappa+0.0)
|
|
error=sqrt(var/Neff)
|
|
return (mean,var,error,kappa)
|
|
|
|
def Block (x, blockfactor):
|
|
N = len(x)
|
|
Nblocks = N / blockfactor
|
|
xb = zeros(Nblocks)
|
|
for i in range(0,Nblocks):
|
|
start = i*blockfactor
|
|
end = (i+1)*blockfactor
|
|
xb[i] = sum(x[start:end])/(blockfactor+0.0)
|
|
return xb
|
|
|
|
def MeanErrorString (mean, error):
|
|
if (mean!=0.0):
|
|
meanDigits = math.floor(math.log(abs(mean))/math.log(10))
|
|
else:
|
|
meanDigits=2
|
|
if (error!=0.0):
|
|
rightDigits = -math.floor(math.log(error)/math.log(10))+1
|
|
else:
|
|
rightDigits=2
|
|
if (rightDigits < 0):
|
|
rightDigits = 0
|
|
formatstr = '%1.' + '%d' % rightDigits + 'f'
|
|
meanstr = formatstr % mean
|
|
errorstr = formatstr % error
|
|
return meanstr + ' +/- ' + errorstr
|
|
# return (meanstr, errorstr)
|
|
|
|
s = loadtxt(sys.argv[1])
|
|
c = s.shape[0]
|
|
if (len(sys.argv) > 3):
|
|
s = s / float(sys.argv[3])
|
|
if len(sys.argv) > 2:
|
|
first = int(sys.argv[2])
|
|
else:
|
|
first = 20
|
|
|
|
E = s[first:c,1];
|
|
Var = s[first:c,2];
|
|
Weight = s[first:c,3];
|
|
Pop = s[first:c,4];
|
|
Etrial = s[first:c,5];
|
|
|
|
|
|
|
|
|
|
N = len(E)
|
|
#print "Blockfactor Std Error"
|
|
#errlist = []
|
|
#for factor in range(10,N/2,10):
|
|
# Eblock = Block(E,factor)
|
|
# E2 = Eblock*Eblock
|
|
# avg = sum (Eblock)/(len(Eblock)+0.0)
|
|
# var = var=sum((Eblock-avg)*(Eblock-avg))/(len(Eblock)+0.0)
|
|
# error = sqrt (var/(len(Eblock)+0.0))
|
|
# errlist.append(error)
|
|
# print " %4d %8.6f" % (factor, error)
|
|
|
|
(avg, var, error, kapp) = Stats(E)
|
|
#error = max(errlist)
|
|
print "E = " + MeanErrorString(avg,error) + (" kappa = %1.3f" % (kapp))
|
|
(avg, var, error, kapp) = Stats(Var)
|
|
print "Var = " + MeanErrorString(avg,error)
|
|
(avg, var, error, kapp) = Stats(Pop)
|
|
print "Pop = " + MeanErrorString(avg,error)
|
|
|
|
|
|
|