diff --git a/pyobs/core/README b/pyobs/core/README new file mode 100644 index 0000000..742bd33 --- /dev/null +++ b/pyobs/core/README @@ -0,0 +1 @@ +uwerr deprecated diff --git a/pyobs/core/libcore.py b/pyobs/core/libcore.py index 6bdd69c..59da65d 100644 --- a/pyobs/core/libcore.py +++ b/pyobs/core/libcore.py @@ -22,9 +22,23 @@ ################################################## +def gamma_fft(Wmax,data1,data2=None): + aux1=numpy.fft.fft(data1) + if (data2 is not None): + if (len(data2)!=len(data1)): + raise + else: + aux2=numpy.fft.fft(data2) + else: + aux2=aux1 + + tmp = numpy.fft.ifft(aux1*numpy.conj(aux2)) + return tmp[0:Wmax+1].real + + def gamma(Wmax,data1,data2=None): N=len(data1) - g = numpy.zeros(Wmax, dtype=numpy.double) + g = numpy.zeros(Wmax+1, dtype=numpy.double) if (data2 is not None): if (len(data2)!=N): @@ -32,13 +46,13 @@ def gamma(Wmax,data1,data2=None): else: data2 = data1 - for t in range(Wmax): + for t in range(Wmax+1): g[t] = data1[0:N-t].dot(data2[t:N]) return g def normalize_gamma(gamma, Wmax, N, R): - n = N-R*numpy.arange(0.,Wmax,1.) + n = N-R*numpy.arange(0.,Wmax+1,1.) nn = 1.0/n return numpy.multiply(gamma, nn) @@ -78,7 +92,7 @@ def find_window(rho, N, Stau, texp=None): flag=0 for W in range(1,Wmax): rho_int = rho_int + rho[W] - tauW = Stau/numpy.log((rho_int+1.)/rho_int) + tauW = Stau/numpy.log(numpy.fabs((rho_int+1.)/rho_int)) gW = numpy.exp(-W/tauW) - tauW/numpy.sqrt(W*N) if (gW<0): Wopt = W diff --git a/pyobs/core/utils.py b/pyobs/core/utils.py index 69edd7a..a064d9d 100644 --- a/pyobs/core/utils.py +++ b/pyobs/core/utils.py @@ -33,7 +33,7 @@ def fill_holes(idx_in, data_in, size=1): N = idx_in[-1] - idx_in[0] + 1 idx = numpy.arange(idx_in[0],idx_in[-1]+1,1) data = numpy.array(numpy.zeros(N*size), dtype=numpy.double) - + norm = float(N)/float(len(data_in)/size) j=0 for i in range(N): @@ -48,36 +48,37 @@ def pad_with_zeros(data, before, after): S = numpy.shape(data) n0 = S[-1] data2 = numpy.zeros( S[0:-1]+(before+n0+after,) ) + norm = float(before+n0+after)/float(n0) for i in range(before, before+n0): data2[:,:,i] = numpy.array(data[:,:,i-before]) - return data2 + return data2 * norm # returns the sorted union of the two lists def union(a,b): - u = set(a) | set(b) - return list(u) + u = sorted(set(a) | set(b)) + return list(u) def double_union(a,b,an,bn): - u = list(set(a) | set(b)) - un = [] - for iu in u: - try: - un.append( an[a.index(iu)] ) - except: - un.append( bn[b.index(iu)] ) - return [u, un] + u = list(sorted(set(a) | set(b))) + un = [] + for iu in u: + try: + un.append( an[a.index(iu)] ) + except: + un.append( bn[b.index(iu)] ) + return [u, un] def contains(a,b): - for i in b: - if (a==b): - return b.index(a) - return -1 + for i in b: + if (a==b): + return b.index(a) + return -1 def piechart(x,l,t): - plt.subplots() - plt.title('Components (%d,%d)' % t) - x2 = numpy.array(x)/numpy.sum(x) - plt.pie(x2,labels=l, autopct='%.0f%%',radius=1.0) - plt.axis('equal') - plt.show() + plt.subplots() + plt.title('Components (%d,%d)' % t) + x2 = numpy.array(x)/numpy.sum(x) + plt.pie(x2,labels=l, autopct='%.0f%%',radius=1.0) + plt.axis('equal') + plt.show() diff --git a/pyobs/core/uwerr.py b/pyobs/core/uwerr.py index 8af0898..deab69f 100644 --- a/pyobs/core/uwerr.py +++ b/pyobs/core/uwerr.py @@ -7,14 +7,17 @@ __all__ = ['uwerr'] -def uwerr(data,ncnfg,plot=False,opts=(),Stau=1.5,W=None): +def uwerr(data,data2,ncnfg,plot=False,opts=(),Stau=1.5,W=None): R=len(ncnfg) N = numpy.sum(ncnfg) Wmax = min(ncnfg)/2 gg = numpy.zeros(Wmax, dtype=numpy.double) for ir in range(R): - gg = gg + gamma(Wmax, data[ir]) + if data2 is None: + gg = gg + gamma_fft(Wmax, data[ir]) + else: + gg = gg + gamma_fft(Wmax, data[ir], data2[ir]) gg2 = normalize_gamma(gg, Wmax, N, R) if (W==None): @@ -26,14 +29,12 @@ def uwerr(data,ncnfg,plot=False,opts=(),Stau=1.5,W=None): gg3 = correct_gamma_bias(gg2, Wopt, N) res = tauint(gg3,Wopt,N) - if (plot==True) or isinstance(plot,str): + if (plot==True): rho = gg3/gg3[0] drho = libcore.compute_drho(rho,N) if (plot==True): - _rho_plotter(rho,drho,Wmax,Wopt,res[1:3],opts) - if isinstance(plot,str): - _rho_plotter(rho,drho,Wmax,Wopt,res[1:3],opts,plot) + _rho_plotter(rho,drho,Wmax,Wopt,res[1:3],opts[:-1],opts[-1]) return res diff --git a/pyobs/ensdata.py b/pyobs/ensdata.py index 62378ff..c2dbb8a 100644 --- a/pyobs/ensdata.py +++ b/pyobs/ensdata.py @@ -1,7 +1,7 @@ import numpy import matplotlib.pyplot as plt -from .vwerr import uwerr, uwerr_texp +from .vwerr import uwerr, uwerr_texp, normac from .core.utils import fill_holes, irregular_measurements from .jerr import jackknife_error @@ -67,11 +67,53 @@ def uwerr(self,plot,pfile,pars=(1.5,None)): dtau = numpy.zeros(self.dims) for i in range(self.dims[0]): for j in range(self.dims[1]): - res = uwerr([rd.data[i,j] for rd in self.rdata], [rd.ncnfg for rd in self.rdata], plot, (self.name,i,j,pfile), pars[0], pars[1]) + res = uwerr([rd.data[i,j] for rd in self.rdata], None, + [rd.ncnfg for rd in self.rdata], + plot, (self.name,i,j,pfile), pars[0], pars[1]) sigma[i,j] = res[0] tau[i,j] = res[1] dtau[i,j] = res[2] return [sigma, tau, dtau] + + + def uwcov(self,plot,pars=(1.5,None)): + if isinstance(plot,str): + raise TypeError('uwcov: saving of plots to file not supported'); + sigma = numpy.zeros(self.dims+self.dims) + tau = numpy.zeros(self.dims+self.dims) + dtau = numpy.zeros(self.dims+self.dims) + for i in range(self.dims[0]): + for k in range(i,self.dims[0]): + for j in range(self.dims[1]): + for l in range(j,self.dims[1]): + res = uwerr([rd.data[i,j] for rd in self.rdata], + [rd.data[k,l] for rd in self.rdata], + [rd.ncnfg for rd in self.rdata], + plot, (self.name,i,j,None), pars[0], pars[1]) + sigma[i,j,k,l] = res[0] + tau[i,j,k,l] = res[1] + dtau[i,j,k,l] = res[2] + if (i!=k) or (j!=l): + sigma[k,l,i,j] = res[0] + tau[k,l,i,j] = res[1] + dtau[k,l,i,j] = res[2] + return [sigma, tau, dtau] + + + def normac(self,Wmax): + rho = numpy.zeros(self.dims+self.dims+(Wmax+1,)) + drho= numpy.zeros(self.dims+self.dims+(Wmax+1,)) + for i in range(self.dims[0]): + for k in range(i,self.dims[0]): + for j in range(self.dims[1]): + for l in range(j,self.dims[1]): + res = normac([rd.data[i,j] for rd in self.rdata], + [rd.data[k,l] for rd in self.rdata], + [rd.ncnfg for rd in self.rdata], Wmax) + rho[i,j,k,l,:] = res[0] + drho[i,j,k,l,:] = res[1] + return [rho,drho] + def uwerr_texp(self,plot,pfile,pars=(1.5,0.0,2,0,None)): sigma = numpy.zeros(self.dims) @@ -89,8 +131,9 @@ def naive_err(self): sigma = numpy.zeros(self.dims) for i in range(self.dims[0]): for j in range(self.dims[1]): - res = uwerr([rd.data[i,j] for rd in self.rdata], [rd.ncnfg for rd in self.rdata], - False, (self.name,i,j,''), 1.5, 0) + res = uwerr([rd.data[i,j] for rd in self.rdata], None, + [rd.ncnfg for rd in self.rdata], + False, (self.name,i,j,None), 1.5, 0) sigma[i,j] = res[0] return sigma @@ -131,21 +174,7 @@ def mean(self,data): def create(self,idx,data,mean): subtraction = numpy.array([mean for _ in range(len(idx))]).flatten() delta = data - subtraction - - # check if data is irregular and fix it - if (irregular_measurements(idx)==True): - tmp = fill_holes(idx, delta) - else: - tmp = [idx, delta] - - self.ncnfg = len(idx) - self.idx = list(tmp[0]) - self.data = numpy.zeros(self.dims+(self.ncnfg,)) - for k in range(self.ncnfg): - for i in range(self.dims[0]): - for j in range(self.dims[1]): - self.data[i,j,k] = tmp[1][ (k*self.dims[1] + i)*self.dims[0] + j] - #self.data = numpy.reshape(tmp[1], (self.ncnfg,)+self.dims[::-1]).T + self.fill(idx,delta) def clone(self,full): res = rdata(self.id,self.name,self.dims) @@ -157,13 +186,26 @@ def clone(self,full): res.data = numpy.zeros(self.dims + (self.ncnfg,)) return res - def fill(self,idx,data=None): - self.ncnfg = len(idx) - self.idx = list(idx) - if (data==None): - self.data = numpy.zeros(self.dims + (self.ncnfg,)) + def fill(self,idx,delta=None): + if delta is None: + # fixes idx if irregular + tmp = [numpy.arange(idx[0],idx[-1]+1,1), delta] else: - self.data = numpy.array(data) + # check if data is irregular and fix it + if (irregular_measurements(idx)==True): + tmp = fill_holes(idx, delta,numpy.prod(self.dims)) + else: + tmp = [idx, delta] + + self.ncnfg = len(tmp[0]) + self.idx = list(tmp[0]) + self.data = numpy.zeros(self.dims+(self.ncnfg,)) + if delta is None: + return + for k in range(self.ncnfg): + for i in range(self.dims[0]): + for j in range(self.dims[1]): + self.data[i,j,k] = tmp[1][ (k*self.dims[0] + i)*self.dims[1] + j] class cdata: diff --git a/pyobs/mathfcts.py b/pyobs/mathfcts.py index aa1f91d..fecd03b 100644 --- a/pyobs/mathfcts.py +++ b/pyobs/mathfcts.py @@ -6,7 +6,7 @@ from sympy import Matrix, diff, lambdify __all__ = ['math_sum','math_det','math_tr', - 'math_inv','math_add','math_mul', + 'math_inv','math_eig','math_add','math_mul', 'math_dot','math_scalar','math_binary'] @@ -102,6 +102,28 @@ def math_inv(mean): return [func, [grad]] +def math_eig(mean): + dims = numpy.shape(mean) + [w,v] = numpy.linalg.eig(mean) + + gradw = numpy.zeros((dims[0],1)+dims) + gradv = numpy.zeros(dims+dims) + dfunc = numpy.zeros(dims) + + for k in range(dims[0]): + for l in range(dims[1]): + dfunc[k,l] = 1.0 + tmp0 = (v.T).dot(dfunc).dot(v) + dfunc[k,l] = 0.0 + for i in range(dims[0]): + gradw[i,0,k,l] = tmp0[i,i] + for j in range(dims[1]): + if i!=j: + gradv[:,i,k,l] += (tmp0[i,j]/(w[i]-w[j]))*v[:,j] + func = [numpy.reshape(w,(dims[0],1)), v] + dfunc = [gradw, gradv] + return [func, dfunc] + ####################################################################### def math_add(mean1,mean2): @@ -138,11 +160,13 @@ def math_dot(mean1,mean2): M2 = SymMat(d2,'M2') s = Subs([mean1,mean2],[M1, M2]) - h = [[0]*d2[1]]*d1[0] + h = [] for i in range(d1[0]): + hh = [0]*d2[1] for j in range(d2[1]): for k in range(d1[1]): - h[i][j] = h[i][j] + M1[i,k]*M2[k,j] + hh[j] = hh[j] + M1[i,k]*M2[k,j] + h.append(hh) expr = Matrix(h) #expr = Matrix(M2.dot(M1)).reshape(d1[0],d2[1]).transpose() diff --git a/pyobs/observa.py b/pyobs/observa.py index fa901bf..feb44f1 100644 --- a/pyobs/observa.py +++ b/pyobs/observa.py @@ -15,7 +15,7 @@ __all__ = ['observa', 'derobs', 'log', 'exp', 'sin', 'cos', 'arcsin', 'arccos', 'sinh', 'cosh', 'arcsinh', 'arccosh', - 'trace', 'det', 'inv'] + 'trace', 'det', 'inv','eig'] class InputError(Exception): pass @@ -1034,6 +1034,32 @@ def inv(x): return derobs([x], f, df) +def eig(x): + """ + Compute the eigenvalues and eigenvectors of a square symmetric matrix + + Parameters + ---------- + x : observa + must be a square matrix + + Returns + ------- + w, v : observa + the eigenvalues and a matrix whose columns + correspond to the eigenvectors + + Examples + -------- + >>> [w, v] = eig(x) + """ + + if (x.dims[0]!=x.dims[1]): + raise InputError('unsupported operation for non-square matrices') + [f, df] = math_eig(x.mean) + return [derobs([x], f[0], [df[0]]), derobs([x], f[1], [df[1]])] + + def derobs(inps,func,dfunc=None): """ Compute the derived observable from a given function. It is highly recommended to avoid the usage of this function @@ -1100,7 +1126,7 @@ def derobs(inps,func,dfunc=None): grad = [] for i in range(len(inps)): if dfunc is None: - grad.append( numerical_derivative(all_mean,func,i,inps[i].naive_err()) ) + grad.append( numerical_derivative(all_mean,func,i,inps[i].naive_err()*1e-4) ) else: grad.append( numpy.array(dfunc[i]) )