fieldpy.core.helper

1 """ 2 A collection of useful helper functions. 3 """ 4 import copy 5 import pdb 6 import cPickle as pickle 7 import subprocess 8 import os 9 10 import numpy as np 11 import matplotlib.mlab as mlab 12 from scipy.interpolate import splprep, splev 13 from scipy.signal import gaussian 14 15

16 -def pickle_(var, filename):

17 """ 18 Pickels var into filename. 19 20 @param var: the varibale to pickle 21 22 @type filename: string 23 @param filename: the name of the picklefile (reccomended to use .pkl extension) 24 """ 25 with open(filename, 'wb') as fil: 26 pickle.dump(var, fil, protocol=pickle.HIGHEST_PROTOCOL)

27

28 -def unpickle(filename):

29 with open(filename, 'rb') as fil: 30 return pickle.load(fil)

31

32 -def nonmasked_values(ma):

33 """ 34 Returns an array containing only the non-masked values of a masked 35 array and the indices of them. 36 """ 37 return np.asarray(ma.data[~ma.mask]), ~st.mask

38

39 -def append_field2struct_arr(struct_ar, name, vec, dtype=None):

40 """ 41 Returns a new structured array with the arrays appended. (the 42 arrays cannot be masked) 43 44 @param struct_ar: structured array to be appended to 45 46 @type name: string 47 @param name: the name of the new datatype 48 49 @type vec: np.ndarray 50 @param vec: the vector to be appended (its length has to be == struct_ar.shape[0]) 51 52 @param dtype: the dtype to be used (defaults to dtype of arr) 53 54 >>> r = np.zeros((3,), dtype=[('foo', int), ('bar', float)]) 55 >>> toap = np.array([1,2,3]) 56 >>> append_field2struct_arr(r, 'new', toap) 57 array([(0, 0.0, 1), (0, 0.0, 2), (0, 0.0, 3)], 58 dtype=[('foo', '<i8'), ('bar', '<f8'), ('new', '<i8')]) 59 """ 60 # check 61 if struct_ar.dtype.names is None: 62 raise(TypeError( 63 "First argument is not a structured array!")) 64 # figure out new dtype 65 if dtype is None: 66 dtype = vec.dtype 67 if dtype.names is not None: 68 raise(TypeError( 69 "vec or dtype is that of a structured array!")) 70 # dtype of new field 71 ndtype = [(name, dtype)] 72 # dtype of whole thing 73 newdtype = np.dtype(struct_ar.dtype.descr + ndtype) 74 new_struct_arr = np.ndarray(struct_ar.shape, dtype=newdtype) 75 for name_ in struct_ar.dtype.names: 76 new_struct_arr[name_] = struct_ar[name_] 77 new_struct_arr[name] = vec 78 return new_struct_arr

79 80

81 -def append_mrec2mrec(mrec, toap_mr):

82 """ 83 Return a new masked structured array with a column appended from 84 another masked structured array. 85 86 @type mrec: record np.ma.core.MaskedArray 87 @param mrec: the masked record array to which will be appended 88 89 @type toap_mr: record np.ma.core.MaskedArray 90 @param toap_mr: record np.ma.core.MaskedArray of the right size to be appended 91 92 @rtype: np.ma.core.MaskedArray 93 @return: The masked record array with appended fields 94 95 @note: Modified from 96 U{http://mail.scipy.org/pipermail/numpy-discussion/2007-September/029357.html} 97 98 >>> r = np.zeros((3,), dtype=[('foo', int), ('bar', float)]) 99 >>> mr = r.view(np.ma.MaskedArray) 100 >>> mr[0] = np.ma.masked 101 >>> toap = np.ones((3,), dtype=[('foobar', int)]) 102 >>> toap_mr = toap.view(np.ma.MaskedArray) 103 >>> toap_mr[2] = np.ma.masked 104 >>> append_mrec2mrec(mr, toap_mr) 105 masked_array(data = [(--, --, 1) (0, 0.0, 1) (0, 0.0, --)], 106 mask = [(True, True, False) (False, False, False) (False, False, True)], 107 fill_value = (999999, 1e+20, 999999), 108 dtype = [('foo', '<i8'), ('bar', '<f8'), ('foobar', '<i8')]) 109 <BLANKLINE> 110 """ 111 # figure out new dtype 112 dtype = toap_mr.dtype 113 newdtype = np.dtype(mrec.dtype.descr + dtype.descr) 114 newdtype_bool = np.dtype([(nt[0], 'bool') for nt in newdtype.descr]) 115 # make a rec_array with the right data in it 116 newcolname = toap_mr.dtype.names[0] 117 data = append_field2struct_arr(mrec, newcolname, toap_mr[newcolname]) 118 # make a new mask 119 newmask = append_field2struct_arr(mrec.mask, newcolname, toap_mr[newcolname].mask) 120 # creat the new np.ma.core.MaskedArray 121 newmrec = np.ma.core.MaskedArray(data, mask=newmask, dtype=newdtype) 122 # for field in mrec.dtype.fields: 123 # newmrec[field] = mrec[field] 124 # newmrec[name] = toap_mr 125 # mrec = newmrec 126 return newmrec

127

128 -def append_ma2mrec(mrec, toap_ma, name):

129 """ 130 Return a new masked record array with a column appended from a 131 masked array. It is assumed that whole records are masked and not 132 individual entires in a record. The mask will the constructed 133 with a locial or operation. 134 135 @type mrec: record np.ma.core.MaskedArray 136 @param mrec: the masked record array to which will be appended 137 138 @type toap_ma: np.ma.core.MaskedArray (no-record) 139 @param toap_ma: np.ma.core.MaskedArray of the right size to be appended 140 141 @type name: string 142 @param name: the name to give to the new colum 143 144 @rtype: np.ma.core.MaskedArray 145 @return: The masked record array with appended fields 146 147 @note: Modified from 148 U{http://mail.scipy.org/pipermail/numpy-discussion/2007-September/029357.html} 149 150 >>> r = np.zeros((3,), dtype=[('foo', int), ('bar', float)]) 151 >>> mr = r.view(np.ma.MaskedArray) 152 >>> mr[0] = np.ma.masked 153 >>> toap = np.ones((3,), dtype=int) 154 >>> toap_ma = toap.view(np.ma.MaskedArray) 155 >>> toap_ma[2] = np.ma.masked 156 >>> append_ma2mrec(mr, toap_ma, 'foobar') 157 masked_array(data = [(--, --, 1) (0, 0.0, 1) (0, 0.0, --)], 158 mask = [(True, True, False) (False, False, False) (False, False, True)], 159 fill_value = (999999, 1e+20, 999999), 160 dtype = [('foo', '<i8'), ('bar', '<f8'), ('foobar', '<i8')]) 161 <BLANKLINE> 162 """ 163 # figure out new dtype 164 dtype = toap_ma.dtype 165 newdtype = np.dtype(mrec.dtype.descr + [(name, dtype.descr[0][1])]) 166 newdtype_bool = np.dtype([(nt[0], 'bool') for nt in newdtype.descr]) 167 # make a rec_array with the right data in it 168 newcolname = name 169 data = append_field2struct_arr(mrec, newcolname, toap_ma) 170 # make a new mask 171 newmask = append_field2struct_arr(mrec.mask, newcolname, toap_ma.mask) 172 # creat the new np.ma.core.MaskedArray 173 newmrec = np.ma.core.MaskedArray(data, mask=newmask, dtype=newdtype) 174 # for field in mrec.dtype.fields: 175 # newmrec[field] = mrec[field] 176 # newmrec[name] = toap_ma 177 # mrec = newmrec 178 return newmrec

179

180 -def append_a2mrec(mrec, toap_a, name):

181 """ 182 Return a new masked record array with a column appended from a 183 array. It is assumed that whole records are masked and not 184 individual entires in a record. The mask will the constructed 185 with a locial or operation. 186 187 @type mrec: record np.ma.core.MaskedArray 188 @param mrec: the masked record array to which will be appended 189 190 @type toap_a: np.array (no-record) 191 @param toap_a: np.array of the right size to be appended 192 193 @type name: string 194 @param name: the name to give to the new colum 195 196 @rtype: np.ma.core.MaskedArray 197 @return: The masked record array with appended fields 198 199 @note: thhs just wraps append_field2struct_arr 200 201 >>> r = np.zeros((3,), dtype=[('foo', int), ('bar', float)]) 202 >>> mr = r.view(np.ma.MaskedArray) 203 >>> mr[0] = np.ma.masked 204 >>> toap = np.ones((3,), dtype=int) 205 >>> append_a2mrec(mr, toap, 'foobar') 206 masked_array(data = [(--, --, 1) (0, 0.0, 1) (0, 0.0, 1)], 207 mask = [(True, True, False) (False, False, False) (False, False, False)], 208 fill_value = (999999, 1e+20, 999999), 209 dtype = [('foo', '<i8'), ('bar', '<f8'), ('foobar', '<i8')]) 210 <BLANKLINE> 211 """ 212 newcolname = name 213 dtype = toap_a.dtype 214 newdtype = np.dtype(mrec.dtype.descr + [(name, dtype.descr[0][1])]) 215 data = append_field2struct_arr(mrec, newcolname, toap_a) 216 newmask = np.zeros(len(toap_a), dtype=bool) 217 newmask = append_field2struct_arr(mrec.mask, newcolname, newmask) 218 newmrec = np.ma.core.MaskedArray(data, mask=newmask, dtype=newdtype) 219 return newmrec

220

221 -def mrec_keep_fields(mrec, names):

222 """ 223 Return a new numpy record array with only fields listed in names 224 225 @note: copied and adjusted from mlab.rec_keep_fields 226 """ 227 data = [] 228 mask = [] 229 dtype = [] 230 for name in names: 231 data.append(mrec[name]) 232 mask.append(mrec.mask[name]) 233 dtype.append((name, mrec[name].dtype)) 234 data = np.rec.fromarrays(data, names=names) 235 mask = np.rec.fromarrays(mask, names=names) 236 dtype = np.dtype(dtype) 237 238 newmrec = np.ma.core.MaskedArray(data, mask=mask, dtype=dtype) 239 return newmrec

240

241 -def filter_gauss(t, x, delta_t, sigma_factor=2.):

242 """Filtes a time signal with a weighted running average in a given 243 time intervall +/-delta_t. 244 The weight is a gaussian such that at t+/-delta_t is at 245 sigma_factor sigma, e.g sigma_factor=2 is c.a. 1/8 of the weight 246 of the center. The filter is symmetric, i.e. no phase shift. 247 248 @param t: time vector 249 @param x: data vector (or array) 250 @param delta_t: filter time interval 251 252 @note: 253 - Really slow as it itterates over all elements... maybe a cython case? 254 - Adapted and refined from my old dataprocessing stuff 255 256 >>> t = np.linspace(0,100,10) 257 >>> x = np.sin(t) 258 >>> filter_gauss(t, x, 10*np.pi) 259 array([ 0. , -0.4581598 , -0.03226275, 0.1324937 , 0.06281049, 260 -0.11801212, -0.09001935, 0.09725727, 0.36815581, -0.50636564]) 261 >>> tm = t.view(np.ma.MaskedArray) 262 >>> xm = x.view(np.ma.MaskedArray) 263 >>> xm[3] = np.ma.masked 264 >>> filter_gauss(tm, xm, 10*np.pi) 265 masked_array(data = [0.0 -0.458159800543 -0.0375228441936 -0.219964900938 0.0730510540235 266 -0.386390056324 -0.090019349563 0.0972572677569 0.368155810038 267 -0.50636564111], 268 mask = [False False False False False False False False False False], 269 fill_value = 1e+20) 270 <BLANKLINE> 271 """ 272 delta_t = float(delta_t) 273 # do a consistency check: 274 if delta_t >= (np.max(t)-np.min(t))/2.: 275 raise(TypeError( 276 "delta_t is bigger than the time interval of t.")) 277 278 if np.ma.is_masked(t): 279 raise ValueError('Masking the time dimension is not supported') 280 # if x is a row vector, make a column vector: 281 row_vec = False 282 if len(x.shape) == 1: 283 x = x.reshape(len(x),1) 284 row_vec = True 285 286 # make a decent kernel 287 # use a gauss shape with width such that the outermost point is at 2 sigma: 288 def gauss(xx, sigma, mu): 289 """gauss function: gauss(x, sigma=1, u=0)""" 290 return np.exp(-(xx-mu)**2/(2.*sigma**2))

291 292 new_x = x.copy() # empty_like is buggy: see numpy trac 293 new_x.fill(0.) 294 if np.ma.is_masked(x): 295 new_x.mask = False 296 # do a loop over all elements 297 for ii in np.arange(0,len(x)): 298 # need some extra stuff for masked arrays 299 # if np.ma.is_masked(t[ii]): 300 # raise ValueError('Masking the time dimension is not supported') 301 # #t_cur = t.data[ii] 302 # else: 303 # 304 t_cur = t[ii] 305 306 # find values which lie in the time range: 307 t_ind = np.where(np.logical_and((t <= t_cur + delta_t), (t >= t_cur - delta_t)))[0] 308 n_ind = len(t_ind) 309 t_ind_u = t_ind[np.where(t_cur < t[t_ind])[0]] 310 t_ind_l = t_ind[np.where(t_cur > t[t_ind])[0]] 311 312 #pdb.set_trace() 313 314 # use values at indices which have a symmetric unmasked sibling: 315 #print '\nind = %i, t_ind = %s' % (ii, str(t_ind)) 316 t_ind_good = np.array([ii]) # dont care whether this is masked as there is no partner 317 for jj in t_ind_l: 318 if np.ma.is_masked(x[jj]): 319 #print '--' 320 continue # do nothing if x[jj] is masked 321 else: 322 t_jj_u = t[t_ind_u] 323 324 # find corresponding t_ind_u, if it exists 325 ind_tmp = t_ind_u[np.where( 326 np.abs(t[jj]-t_cur - (-t_jj_u+t_cur)) < delta_t/n_ind )[0]] 327 if np.ma.is_masked(x[ind_tmp]): 328 #print '++' 329 continue # dicharge if masked 330 # print ind_tmp 331 332 if ind_tmp.shape[0] == 1: 333 #print 'in', ind_tmp, t_ind, t_ind_l 334 t_ind_good = np.concatenate((ind_tmp,t_ind_good)) 335 t_ind_good = np.concatenate((np.array([jj]),t_ind_good)) 336 t_ind_good.sort() 337 # print t_ind_good 338 339 weights = gauss(t[t_ind_good], delta_t/sigma_factor, t_cur) 340 if np.ma.is_masked(x[t_ind_good]): 341 weights[np.where(x[t_ind_good].mask)[0]] = np.ma.masked 342 # print weights 343 344 # if x.mask[ii]==False: 345 # pdb.set_trace() 346 #pdb.set_trace() 347 348 # normalize to one: 349 weights = weights/weights.sum() 350 # now calcutate value: 351 new_x[ii] = (weights.reshape(weights.shape[0],1)* x[t_ind_good,...]).sum(axis=0) 352 353 if row_vec: 354 new_x = new_x.ravel() 355 return new_x 356

357 -def filter_gauss_conv(t, x, delta_t, sigma_factor=2.0):

358 """ Filters a regularly-sampled signal with a Gaussian-weighted 359 average with half-window length delta_t. The Gaussian bell is 360 chosen such that at +/-delta_t it is valued at sigma_factor, e.g. 361 sigma_factor=2 is c.a. 1/8 of the weight of the center. There 362 is no phase shift. 363 364 Warning: If x is masked, then NaN-values will be laundered into 365 masked values in the output. This means that if x contains both a 366 mask and NaNs, the output will only have a mask, extended over the 367 NaNs. 368 369 @param t: time vector 370 @param x: data vector (or array) 371 @param delta_t: filter time interval 372 373 @note: 374 - Performs filtering using a np.convolve. 375 376 >>> t = np.linspace(0,100,10) 377 >>> x = np.sin(t) 378 >>> filter_gauss_conv(t, x, 10*np.pi) 379 array([-0.35491283, -0.25778845, -0.06943494, 0.07441436, 0.03527717, 380 -0.06628086, -0.05055887, 0.10013494, 0.20257496, 0.20875278]) 381 >>> tm = t.view(np.ma.MaskedArray) 382 >>> xm = x.view(np.ma.MaskedArray) 383 >>> xm[3] = np.ma.masked 384 >>> filter_gauss_conv(tm, xm, 10*np.pi) 385 masked_array(data = [-- -- -- 0.0380257449674 0.035277166102 -0.066280860855 -0.0505588746925 386 0.100134936194 0.202574957573 0.208752775728], 387 mask = [ True True True False False False False False False False], 388 fill_value = 1e+20) 389 <BLANKLINE> 390 """ 391 delta_t = float(delta_t) 392 393 # Make sure that dt < range(t) 394 if delta_t >= (np.max(t)-np.min(t)) / 2.0: 395 # I changed this to a ValueError, but it could go back to 396 # TypeError if it's caught somewhere else -njw 397 raise(ValueError( 398 "delta_t is must be smaller than half the range of t")) 399 if np.ma.is_masked(t): 400 raise TypeError("Masking the time dimension is not supported") 401 402 # Relate the sampling interval to dt 403 # Use a mean just in case t isn't strictly regularly spaced 404 halfwin = int(round( delta_t / np.diff(t).mean() )) 405 406 # Get a kernel - just grap one from scipy 407 kernel = gaussian(2*halfwin+1, halfwin/sigma_factor) 408 kernel /= kernel.sum() 409 410 # Deal with masked x 411 if np.ma.is_masked(x): 412 was_masked = True 413 x_nan = x.data 414 x_nan[x.mask is True] = np.nan 415 else: 416 was_masked = False 417 x_nan = x 418 419 # Interpolate over NaN values 420 if True in np.isnan(x_nan): 421 x_valid = x_nan[np.nonzero(np.isnan(x_nan)==False)] 422 t_valid = t[np.nonzero(np.isnan(x_nan)==False)] 423 x_interp = np.interp(t, t_valid, x_valid) 424 else: 425 x_interp = x_nan 426 427 # Pad the data array 428 pad0 = x[:halfwin].mean() 429 padf = x[-halfwin:].mean() 430 x_padded = np.hstack([pad0*np.ones(halfwin), x_interp, padf*np.ones(halfwin)]) 431 432 # Perform the filtering 433 x_filtered = np.convolve(x_padded, kernel, mode='valid') 434 435 if was_masked: 436 # Build a new masked array around x_filtered in-place 437 x_filtered = np.ma.masked_where(np.isnan(x_filtered), 438 x_filtered) 439 440 return x_filtered

441

442 -def filter_and_resample_slow(time, data, filter_window, 443 new_sample_int, first_sample=None):

444 """Filters and resamples a time signal with a running average in a 445 given time intervall +/-filter_window/2. The filter has a PHASE 446 SHIFT for a UN-equally sampled signal. Accepts masked arrays but 447 not "structured arrays" or "recarrays". 448 449 If you have cython (or just build the modul) use the much faster 450 filter_and_resample in helper_fns_cython. 451 452 Having said this, this python implementation is probably really slow! 453 454 @param time: the times of the timeseries 455 @param data: the data of the timeseries 456 @param filter_window: filter time interval 457 @param new_sample_int: the new sampling interval (in units of time) 458 @param first_sample: the time of the first sample in the resampled 459 signal. If None (default) it is 460 (approximatly) set to time[0]+filter_window 461 462 @rtype: list of np.array 463 @return: Returns [new_time, new_data] 464 465 @note: 1) ignores any masks in the time variable. 466 2) will not work if it finds no values within filter window 467 @todo: Correct above point 2) (as well in helper_fns_cython.pyx) 468 >>> tim = np.linspace(-3,90,1000); data = tim**2 469 >>> filter_window = 30; new_sample_int = 30 470 >>> filter_and_resample_slow(tim, data, filter_window, new_sample_int) 471 (array([ 0., 30., 60., 90.]), array([ 62.98183318, 975.90237234, 3678.22848073, 6826.19355542])) 472 >>> data = np.ma.MaskedArray(data) 473 >>> data[1:10] = np.ma.masked 474 >>> filter_and_resample_slow(tim, data, filter_window, new_sample_int) 475 (array([ 0., 30., 60., 90.]), array([ 65.73049119, 975.90237234, 3678.22848073, 6826.19355542])) 476 """ 477 filter_window = filter_window/2. 478 if first_sample is None: 479 first_sample = time[0] + filter_window 480 # make it such that it lies nicely with the sampling interval 481 first_sample = (np.round(first_sample/new_sample_int)) * new_sample_int 482 new_t = np.arange( first_sample, time[-1]+1e-10*time[-1], new_sample_int) 483 new_d = np.empty_like(new_t) 484 ind1 = 0 485 ind2 = 0 486 for ii in range(len(new_t)): 487 ind1 = get_ind_closest(time, new_t[ii]-filter_window, ind1) 488 ind2 = get_ind_closest(time, new_t[ii]+filter_window, ind2) 489 if ind1==-1 or ind2==-1: 490 raise(ValueError('Index not found')) 491 new_d[ii] = data[ind1:ind2+1].mean() 492 493 return new_t, new_d

494 495 496

497 -def get_ind_closest(vec, value, guess_ind=None):

498 """ 499 Finds the index (ind) of vec such that vec[ind]-value is smallest. Only 500 works reliably for monotone vectors. As start value it uses the 501 value found in ind[0]. 502 503 This function is somewhat slower than get_ind_closest but can be 504 called directly from Python. (When calling from cython use 505 get_ind_closest_faster) 506 507 @param vec: the vector, usually a time 508 @param value: the value which is compared to vec 509 @param guess_ind: a guess for index around which the search 510 (symmetric) will be started. 511 512 @return: The number of indices away from 513 514 @note: The algorithm is fairly dum, something like bisection might 515 work much faster. 516 517 >>> tim =np.linspace(-3,100,1000) 518 >>> get_ind_closest(tim, 49.2) 519 506 520 >>> get_ind_closest(tim, 49.2, 400) 521 506 522 """ 523 # checks 524 len_v = len(vec) 525 if guess_ind < 0: 526 guess_ind = 0 527 if guess_ind>len_v-1: 528 guess_ind = len_v-1 529 # figure out sign of the difference between vec[guess_ind] and value 530 si = np.int(np.sign(vec[guess_ind] - value)) 531 532 # do a while loop until the sign changes 533 ii = 0 534 while True: 535 if guess_ind+ii<len_v: 536 s1 = (np.sign(vec[guess_ind+ii] - value)) 537 if guess_ind-ii>=0: 538 s2 = (np.sign(vec[guess_ind-ii] - value)) 539 if s1!=si: 540 # break if sign changes 541 ind = guess_ind + ii 542 d1 = np.abs(vec[ind] - value) 543 d2 = np.abs(vec[ind-1] - value) 544 if d1<d2: 545 return ind 546 else: 547 return ind-1 548 if s2!=si: 549 # break if sign changes 550 ind = guess_ind - ii 551 d1 = vec[ind] - value 552 d2 = vec[ind+1] - value 553 if d1<d2: 554 return ind 555 else: 556 return ind+1 557 if guess_ind+ii>=len_v and guess_ind-ii<0: 558 if abs(guess_ind-len_v) < abs(guess_ind): 559 return len_v-1 560 else: 561 return 0 562 ii += 1

563

564 -def spline_interpolation(vecs, t_par, smoothness=None, spline_order=3, 565 n_knots=None, n_knots_estimate=-1):

566 """Finds a B-spline representation of an N-dimensional curve. Uses 567 scipy.interpolate.fitpack.splprep. 568 569 Given a list of N rank-1 arrays, vecs, which represent a curve in 570 N-dimensional space parametrized by t_par, find a smooth 571 approximating spline curve g(t_par). Uses the FORTRAN routine 572 parcur from FITPACK. 573 574 Example: 575 t = np.hstack((np.linspace(0,4*np.pi, 2000), np.linspace(4*np.pi+0.1, 8*np.pi, 500))) 576 xx = np.sin(t) 577 x = xx + np.random.normal(scale=0.1, size=t.shape) 578 fn,tckp = spline_interpolation([x],t, smoothness=40.0) 579 plot(t,x),hold(True), plot(t, fn(t)[0],'r'), plot(t,xx, 'g') 580 581 @param vecs: A list of sample vector arrays representing the curve. 582 @param t_par: An array of parameter values (probably time). 583 @param smoothness: smoothness parameter (default determinded by splprep) 584 @parameter spline_order: spline order (default & reccomended 3) 585 586 @note: adapted from 587 http://www.scipy.org/Cookbook/Interpolation#head-34818696f8d7066bb3188495567dd776a451cf11 588 """ 589 # find the knot points 590 tckp,u = splprep(vecs, u=t_par, s=smoothness, k=spline_order, 591 nest=n_knots_estimate) 592 593 # return a function which evaluates the spline 594 return lambda tt : splev(tt, tckp), tckp

595 596 597

598 -def polyval(p, x):

599 """ 600 Copied from numpy and adjusted to work for masked arrays. 601 602 Evaluate a polynomial at specific values. 603 604 If `p` is of length N, this function returns the value: 605 606 ``p[0]*x**(N-1) + p[1]*x**(N-2) + ... + p[N-2]*x + p[N-1]`` 607 608 If `x` is a sequence, then `p(x)` is returned for each element of `x`. 609 If `x` is another polynomial then the composite polynomial `p(x(t))` 610 is returned. 611 612 Parameters 613 ---------- 614 p : array_like or poly1d object 615 1D array of polynomial coefficients (including coefficients equal 616 to zero) from highest degree to the constant term, or an 617 instance of poly1d. 618 x : array_like or poly1d object 619 A number, a 1D array of numbers, or an instance of poly1d, "at" 620 which to evaluate `p`. 621 622 Returns 623 ------- 624 values : ndarray or poly1d 625 If `x` is a poly1d instance, the result is the composition of the two 626 polynomials, i.e., `x` is "substituted" in `p` and the simplified 627 result is returned. In addition, the type of `x` - array_like or 628 poly1d - governs the type of the output: `x` array_like => `values` 629 array_like, `x` a poly1d object => `values` is also. 630 631 See Also 632 -------- 633 poly1d: A polynomial class. 634 635 Notes 636 ----- 637 Horner's scheme [1]_ is used to evaluate the polynomial. Even so, 638 for polynomials of high degree the values may be inaccurate due to 639 rounding errors. Use carefully. 640 641 References 642 ---------- 643 .. [1] I. N. Bronshtein, K. A. Semendyayev, and K. A. Hirsch (Eng. 644 trans. Ed.), *Handbook of Mathematics*, New York, Van Nostrand 645 Reinhold Co., 1985, pg. 720. 646 647 Examples 648 -------- 649 >>> np.polyval([3,0,1], 5) # 3 * 5**2 + 0 * 5**1 + 1 650 76 651 >>> np.polyval([3,0,1], np.poly1d(5)) 652 poly1d([ 76.]) 653 >>> np.polyval(np.poly1d([3,0,1]), 5) 654 76 655 >>> np.polyval(np.poly1d([3,0,1]), np.poly1d(5)) 656 poly1d([ 76.]) 657 658 """ 659 p = np.asarray(p) 660 if isinstance(x, np.poly1d): 661 y = 0 662 else: 663 y = np.zeros_like(x) 664 for i in range(len(p)): 665 y = x * y + p[i] 666 return y

667

668 -def get_current_fieldpy_git_hash():

669 """ 670 Returns the has of the current git revision of the fieldpy 671 package. 672 673 @return: Git SHA1 hash 674 """ 675 dir_of_this_file = os.path.dirname(__file__) 676 return get_current_git_hash(dir_of_this_file)

677

678 -def get_current_git_hash(dirname):

679 """ 680 Returns the has of the current git revision of directory passed in. 681 682 @param dirname: directory in which the git revision is queried 683 684 @return: Git SHA1 hash 685 """ 686 git_cmd = 'git rev-parse --verify HEAD' 687 cmd_str = 'cd ' + dirname +';' + git_cmd # cd into directory of this file 688 try: 689 # return subprocess.check_output(cmd_str, shell=True).strip() 690 return _check_output(cmd_str, shell=True).strip() 691 except subprocess.CalledProcessError as err: 692 return 'git version control not running or producing error: no current version'

693 # except: 694 # return 'other error processing git' 695 696

697 -def _check_output(*popenargs, **kwargs):

698 r""" 699 Copied from python 2.7 standard libary as python 2.6 doesn't have 700 it. Remove after Flavien updated to 2.7. 701 702 Run command with arguments and return its output as a byte string. 703 704 If the exit code was non-zero it raises a CalledProcessError. The 705 CalledProcessError object will have the return code in the returncode 706 attribute and output in the output attribute. 707 708 The arguments are the same as for the Popen constructor. Example: 709 710 >>> _check_output(["ls", "-l", "/dev/null"]) # doctest:+ELLIPSIS 711 'crw-rw-rw- 1 root root 1, ... 712 713 The stdout argument is not allowed as it is used internally. 714 To capture standard error in the result, use stderr=STDOUT. 715 716 >>> _check_output(["/bin/sh", "-c", 717 ... "ls -l non_existent_file ; exit 0"], 718 ... stderr=subprocess.STDOUT) 719 'ls: cannot access non_existent_file: No such file or directory\n' 720 """ 721 if 'stdout' in kwargs: 722 raise ValueError('stdout argument not allowed, it will be overridden.') 723 process = subprocess.Popen(stdout=subprocess.PIPE, *popenargs, **kwargs) 724 output, unused_err = process.communicate() 725 retcode = process.poll() 726 if retcode: 727 cmd = kwargs.get("args") 728 if cmd is None: 729 cmd = popenargs[0] 730 raise subprocess.CalledProcessError(retcode, cmd, output=output) 731 return output

732

Source Code for Module fieldpy.core.helper_fns