snobfit.util

1 """ 2 Supplementary or Auxiliary Functions for Snobfit 3 4 Most of them are the functions of the Rough MATLAB-Numpy Equivalents. 5 """ 6 7 import numpy 8 import scipy 9 import scipy.linalg 10 11 12 # Some constants 13 eps = numpy.finfo('d').eps 14 inf = numpy.inf 15 nan = numpy.nan 16 17 # Some internal functions of numpy 18 rand = numpy.random.rand 19 norm = numpy.linalg.norm 20 inv = numpy.linalg.inv 21 triu = numpy.triu 22 isnan = numpy.isnan 23 24 25 # ----------------------------------------------------------

26 -def within(low, x, high):

27 return numpy.logical_and(low <= x, high >= x)

28 29

30 -def feval(funcName, *args):

31 return eval(funcName)(*args)

32 33

34 -def find(*args,**kw):

35 tup = numpy.where(*args,**kw) 36 if len(tup)==0: 37 return numpy.array([],'int') 38 else: 39 return tup[0]

40 41

42 -def toCol( v ):

43 """ 44 v is a vector 45 """ 46 return numpy.matrix( v.reshape(v.size,1) )

47 48

49 -def toRow( m ):

50 """ 51 m is a n x 1 matrix 52 """ 53 n = m.shape[0] 54 return numpy.array( numpy.asarray(m).reshape(1,n)[0] )

55 # numpy.array( m.A.reshape(1,n)[0] ) 56 57

58 -def mldiv(a,b):

59 """ 60 Backslash or left matrix divide. 61 a\b is the matrix division of A into B, which is roughly the 62 same as INV(A)*B 63 """ 64 if a.shape[0] == a.shape[1]: 65 return scipy.linalg.solve(a,b) 66 else: 67 return numpy.asarray( scipy.linalg.lstsq(a,b)[0] ).T[0]

68 69

70 -def chol( a ):

71 """ 72 a is a mat 73 """ 74 n = a.shape[0] 75 return scipy.linalg.cholesky(a)

76 77

78 -def svd( a ):

79 """ 80 a is a mat 81 """ 82 (U,S,V) = scipy.linalg.svd(a,0) 83 return (U,S,V.T)

84 85

86 -def qr(a):

87 (Q,R) = numpy.linalg.qr(a) 88 return R

89 90

91 -def vector(n):

92 """ 93 Return a vector of the given length. 94 """ 95 return numpy.empty(n,'d')

96

97 -def ivector(n=0):

98 """ 99 Return a int vector of the given length. 100 """ 101 return numpy.zeros(n,'int')

102 103

104 -def sort(x):

105 ind = numpy.argsort(x); 106 return (x[ind], ind)

107 108

109 -def std(x):

110 """ 111 STD(X) normalizes by (N-1) where N is the sequence length. 112 This makes STD(X).^2 the best unbiased estimate of the variance 113 if X is a sample from a normal distribution. 114 """ 115 return numpy.std(x,ddof=1)

116 117

118 -def max_(x):

119 """ 120 [Y,I] = MAX(X) returns the indices of the maximum values in vector I. 121 If the values along the first non-singleton dimension contain more 122 than one maximal element, the index of the first one is returned. 123 """ 124 if len(x.shape) > 1: 125 print "It's a Matrix" 126 127 idx = numpy.argmax(x) 128 return ( x[idx], idx )

129 130

131 -def min_(x):

132 """ 133 [Y,I] = MIN(X) returns the index of the minimum value in vector I. 134 If the values along the first non-singleton dimension contain more 135 than one minimal element, the index of the first one is returned. 136 """ 137 if len(x.shape) > 1: 138 print "It's a Matrix" 139 idx = numpy.argmin(x) 140 return ( x[idx], idx )

141 142

143 -def toInt(x):

144 # Make sure it is vector 145 if x.shape[0]==1: 146 x = numpy.reshape( numpy.asarray(x), x.shape[1] ) 147 148 Ix = numpy.array( x, 'int' ) 149 # 10000 is big enough??? Fix me later. 150 151 Ix[ numpy.where(Ix<0) ] = 10000 152 153 return Ix

154 155

156 -def isEmpty(x):

157 if len(x)==0: return True 158 if x.shape[0]==0: return True 159 if x.shape[1]==0: return True 160 return False

161 162

163 -def removeByInd(a, ind):

164 n = len(a) - len(ind) 165 ret = numpy.zeros( n, 'int' ) 166 k = 0 167 for i in xrange(len(a)): 168 if (i not in ind): 169 ret[k] = a[i] 170 k += 1 171 return ret

172 173

174 -def duplicate(v, n, axis=0):

175 #return numpy matrix 176 return numpy.mat( v ).repeat(n, axis=axis)

177 178

179 -def dup(v, n, axis=0):

180 # Return a numpy array 181 if axis == 0: 182 return numpy.array( [v] ).repeat(n, axis=axis) 183 else: 184 return numpy.array( v ).repeat(n, axis=axis)

185 186

187 -def dot3(a,b,c):

188 # Calculate a*(b*c) 189 return numpy.dot(a, numpy.dot(b,c) )

190 191

192 -def crdot(m, v):

193 n = len(v) 194 _sum = 0.0 195 for i in xrange(n): 196 _sum += m[i,0]*v[i] 197 198 return _sum

199 200 #---------------------------------------------------------------

201 -def rsort(x,w=None):

202 """ 203 Sort x in increasing order, remove multiple entries, 204 and adapt weights w accordingly x and w must both be a row or a column 205 default input weights are w=1 206 207 If w==None, the weighted empirical cdf is computed at x 208 dof = len(x) at input is the original number of degrees of freedom 209 210 Warning: when you use this function, make sure x and w is row vector 211 """ 212 if w is None: 213 w = numpy.ones( x.shape ) 214 215 ind = numpy.argsort(x) 216 x = x[ind] 217 w = w[ind] 218 219 # Remove repetitions 220 n = len(x) 221 222 xnew = numpy.append( x[1:n], inf ) 223 224 ind = find( xnew != x ) 225 nn = len(ind) 226 227 x=x[ind] 228 229 ww = numpy.zeros(nn) 230 ww[0] = sum( w[0:ind[0]+1 ] ) 231 for i in xrange( nn-1 ): 232 ww[i+1] = sum( w[ ind[i]:ind[i+1]] ) 233 234 # Get cumulative sum of weights 235 cdfx = numpy.cumsum(ww) 236 237 # Adjust for jumps and normalize 238 cdfx = (cdfx-0.5*ww)/cdfx[nn-1] 239 dof = n 240 241 return (x,ww,cdfx,dof)

242 243 244 245 # --------------------------------------------------------------- 246 if __name__ == '__main__': 247 print rand(5) 248

Source Code for Module snobfit.util