added eig for symmetric matrix; fixed bug in dot

pyobs/core/README

@@ -0,0 +1 @@
+uwerr deprecated

pyobs/core/libcore.py

@@ -22,23 +22,37 @@
 ##################################################
 
 
+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):
             raise
     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

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()

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
 

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:

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()
 

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]) )