1#!/usr/local/bin/python3.8 2 3# python3 status: started 4 5# afni_util.py : general utilities for python programs 6 7# ------------------------------------------------------ 8# no longer usable as a main: see afni_python_wrapper.py 9# ------------------------------------------------------ 10 11import sys, os, math, copy 12from afnipy import afni_base as BASE 13from afnipy import lib_textdata as TD 14import glob 15import pdb 16import re 17 18# global lists for basis functions 19basis_known_resp_l = ['GAM', 'BLOCK', 'dmBLOCK', 'dmUBLOCK', 'SPMG1', 20 'WAV', 'MION'] 21basis_one_regr_l = basis_known_resp_l[:] 22basis_one_regr_l.append('MION') 23stim_types_one_reg = ['file', 'AM1', 'times'] 24g_valid_shells = ['csh','tcsh','sh','bash','zsh'] 25 26# this file contains various afni utilities 17 Nov 2006 [rickr] 27 28def change_path_basename(orig, prefix='', suffix='', append=0): 29 """given a path (leading directory or not) swap the trailing 30 filename with the passed prefix and suffix 31 e.g. C_P_B('my/dir/pickles.yummy','toast','.1D') 32 --> 'my/dir/toast.1D' 33 or with append... 34 e.g. C_P_B('my/dir/pickles.yummy','toast','.1D', append=1) 35 --> 'my/dir/toastpickles.yummy.1D' 36 or maybe we want another dot then 37 e.g. C_P_B('my/dir/pickles.yummy','toast.','.1D', append=1) 38 --> 'my/dir/toast.pickles.yummy.1D' 39 """ 40 if not orig: return '' 41 if not prefix and not suffix: return orig 42 43 (head, tail) = os.path.split(orig) 44 if append: tail = '%s%s%s' % (prefix, tail, suffix) 45 else: tail = '%s%s' % (prefix, suffix) 46 47 if head == '': return tail 48 return "%s/%s" % (head, tail) 49 50# write text to a file 51def write_text_to_file(fname, tdata, mode='w', wrap=0, wrapstr='\\\n', exe=0): 52 """write the given text (tdata) to the given file 53 fname : file name to write (or append) to 54 dtata : text to write 55 mode : optional write mode 'w' or 'a' [default='w'] 56 wrap : optional wrap flag [default=0] 57 wrapstr : optional wrap string: if wrap, apply this string 58 exe : whether to make file executable 59 60 return 0 on success, 1 on error 61 """ 62 63 if not tdata or not fname: 64 print("** WTTF: missing text or filename") 65 return 1 66 67 if wrap: tdata = add_line_wrappers(tdata, wrapstr) 68 69 if fname == 'stdout' or fname == '-': 70 fp = sys.stdout 71 elif fname == 'stderr': 72 fp = sys.stderr 73 else: 74 try: 75 fp = open(fname, mode) 76 except: 77 print("** failed to open text file '%s' for writing" % fname) 78 return 1 79 80 fp.write(tdata) 81 82 if fp != sys.stdout and fp != sys.stderr: 83 fp.close() 84 if exe: 85 try: code = eval('0o755') 86 except: code = eval('0755') 87 try: 88 os.chmod(fname, code) 89 except: 90 print("** failed chmod 755 on %s" % fname) 91 92 return 0 93 94def wrap_file_text(infile='stdin', outfile='stdout'): 95 """make a new file with line wrappers 14 Mar 2014 96 97 The default parameters makes it easy to process as a stream: 98 99 cat INPUT | afni_python_wrapper.py -eval 'wrap_file_text()' 100 or 101 afni_python_wrapper.py -eval 'wrap_file_text()' < INPUT > OUTPUT 102 or 103 afni_python_wrapper.py -eval 'wrap_file_text("INPUT", "OUTPUT")' 104 or 105 afni_python_wrapper.py -eval "wrap_file_text('$f1', '$f2')" 106 """ 107 108 tdata = read_text_file(fname=infile, lines=0, strip=0) 109 if tdata != '': write_text_to_file(outfile, tdata, wrap=1) 110 111 112def read_text_file(fname='stdin', lines=1, strip=1, noblank=0, verb=1): 113 """return the text text from the given file as either one string 114 or as an array of lines""" 115 116 if fname == 'stdin' or fname == '-': fp = sys.stdin 117 else: 118 try: fp = open(fname, 'r') 119 except: 120 if verb: print("** read_text_file: failed to open '%s'" % fname) 121 if lines: return [] 122 else: return '' 123 124 if lines: 125 tdata = fp.readlines() 126 if strip: tdata = [td.strip() for td in tdata] 127 if noblank: tdata = [td for td in tdata if td != ''] 128 else: 129 tdata = fp.read() 130 if strip: tdata.strip() 131 132 fp.close() 133 134 return tdata 135 136def read_top_lines(fname='stdin', nlines=1, strip=0, verb=1): 137 """use read_text_file, but return only the first 'nlines' lines""" 138 tdata = read_text_file(fname, strip=strip, verb=verb) 139 if nlines != 0: tdata = tdata[0:nlines] 140 return tdata 141 142def read_text_dictionary(fname, verb=1, mjdiv=None, mndiv=None, compact=0, 143 qstrip=0): 144 """this is the same as read_text_dict_list(), but it returns a dictionary 145 146 if compact, collapse single entry lists 147 if qstrip, strip any containing quotes 148 """ 149 rv, ttable = read_text_dict_list(fname, verb=verb, mjdiv=mjdiv, mndiv=mndiv, 150 compact=compact, qstrip=qstrip) 151 if rv: return rv, {} 152 153 rdict = {} 154 for row in ttable: 155 if len(row) != 2: 156 print("** RT_dict: table '%s' has bad row length %d" \ 157 % (fname, len(row))) 158 return 1, {} 159 rdict[row[0]] = row[1] 160 161 return 0, rdict 162 163def read_text_dict_list(fname, verb=1, mjdiv=None, mndiv=None, compact=0, 164 qstrip=0): 165 """read file as if in a plain dictionary format (e.g. LABEL : VAL VAL ...) 166 167 mjdiv : major divider can be a single string (':') or a list of them 168 to try, in order 169 (if None, partitioning will be attempted over a few cases) 170 mkdiv : minor divider should be a single string ('SPACE', ',,') 171 172 if either divider is 'SPACE', the natural .split() will be used 173 174 compact: collapse any single entry lists to return 175 DLIST enties of form [LABEL, VAL] or [LABEL, [V0, V1...]] 176 177 qstrip: strip any surrounding quotes 178 179 return status, DLIST 180 where status == 0 on success 181 DLIST is a list of [LABEL, [VAL,VAL...]] 182 183 The list is to keep the original order. 184 """ 185 tlines = read_text_file(fname, lines=1, strip=1, noblank=1, verb=verb) 186 nrows = len(tlines) 187 if nrows == 0: return 0, [] 188 189 # note the major divider 190 if mjdiv == None: 191 major_list = [':', '=', 'SPACE', '::'] 192 elif type(mjdiv) == list: 193 major_list = mjdiv 194 else: 195 major_list = [mjdiv] 196 197 # see if we have a major partitioning 198 for mdiv in major_list: 199 ttable = [] # keep the sorting 200 good = 1 201 for line in tlines: 202 if mdiv == 'SPACE': entries = line.split() 203 else: entries = line.split(mdiv) 204 if len(entries) == 2: 205 entries = [e.strip() for e in entries] 206 ttable.append(entries) 207 else: 208 good = 0 209 break 210 if good: 211 if verb > 1: print("text_dict_list: using major delim '%s'" % mdiv) 212 break 213 214 # did we find a major partitioning? 215 if not good: 216 if verb: 217 print("** failed read_text_dictionary('%s')" % fname) 218 return 1, [] 219 220 # try to guess the minor divider 221 if mndiv == None: 222 has_cc = 0 223 has_c = 0 224 has_s = 0 225 for tline in ttable: 226 if ',,' in tline: 227 has_cc += 1 228 mndiv = ',,' 229 # give double comma priority, since it should be uncommon 230 break 231 if ',' in tline: 232 has_c += 1 233 if mndiv == None: 234 if has_cc: mndiv = ',' 235 else: mndiv = 'SPACE' 236 237 # now actually make a dictionary 238 outtable = [] 239 for tline in ttable: 240 if qstrip: 241 tline[1] = tline[1].strip("'") 242 tline[1] = tline[1].strip('"') 243 if mndiv == 'SPACE': entries = tline[1].split() 244 else: entries = tline[1].split(mndiv) 245 if compact and len(entries) == 1: 246 entries = entries[0] 247 outtable.append([tline[0], entries]) 248 249 return 0, outtable 250 251def convert_table2dict(dlist): 252 """try to convert a table to a dictionary 253 each list entry must have exactly 2 elements 254 255 return status, dictionary 256 """ 257 if type(dlist) == dict: 258 return 0, dlist 259 if not type(dlist) == list: 260 print("** convert_table2dict: unknown input of type %s" % type(dlist)) 261 return 1, None 262 263 rdict = {} 264 for dentry in dlist: 265 if not type(dentry) == list: 266 print("** convert_table2dict: unknown entry type %s" % type(dentry)) 267 return 1, None 268 if len(dentry) != 2: 269 print("** convert_table2dict: invalid entry length %d" % len(dentry)) 270 return 1, None 271 272 # finally, the useful instruction 273 rdict[dentry[0]] = dentry[1] 274 275 return 0, rdict 276 277def write_data_as_json(data, fname='stdout', indent=3, sort=1, newline=1, 278 table2dict=0): 279 """dump to json file; check for stdout or stderr 280 table2dict: if set, convert table Nx2 table into a dictionary 281 return 0 on success 282 """ 283 # possibly convert data to a dictionary 284 if table2dict: 285 rv, dtable = convert_table2dict(data) 286 if rv: 287 print("** write_data_as_json: failed to create file '%s'" % fname) 288 return 1 289 data = dtable 290 291 # import locally, unless it is needed in at least a few functions 292 # (will make default in afni_proc.py, so be safe) 293 try: 294 import json 295 except: 296 print("** afni_util.py: 'json' python library is missing") 297 return 1 298 299 if fname == 'stdout' or fname == '-': 300 fp = sys.stdout 301 elif fname == 'stderr': 302 fp = sys.stderr 303 else: 304 try: fp = open(fname, 'w') 305 except: 306 print("** write_as_json: could not open '%s' for writing" % fname) 307 return 1 308 309 # actual write 310 json.dump(data, fp, sort_keys=sort, indent=indent) 311 312 if newline: fp.write('\n') 313 if fp != sys.stdout and fp != sys.stderr: 314 fp.close() 315 316 return 0 317 318def read_json_file(fname='stdin'): 319 try: import json 320 except: return {} 321 textdata = read_text_file(fname, lines=0) 322 return json.loads(textdata) 323 324def print_dict(pdict, fields=[], nindent=0, jstr=' ', verb=1): 325 istr = ' '*nindent 326 nfields = len(fields) 327 for kk in pdict.keys(): 328 if nfields > 0 and kk not in fields: 329 continue 330 val = pdict[kk] 331 if type(val) == dict: 332 if verb: print('%s%s :' % (istr, kk)) 333 # recur 334 print_dict(val, nindent=(nindent+3), verb=verb) 335 elif type(val) == list: 336 lstr = jstr.join('%s' % v for v in val) 337 if verb: print('%s%s : %s' % (istr, kk, lstr)) 338 else: print(lstr) 339 else: 340 if verb: print('%s%s : %s' % (istr, kk, val)) 341 else: print(val) 342 343def print_json_dict(fname='stdin', fields=[], verb=1): 344 jdata = read_json_file(fname) 345 print_dict(jdata, fields=fields, verb=verb) 346 347def read_AFNI_version_file(vdir='', vfile='AFNI_version.txt', delim=', '): 348 """read AFNI_version.txt from vdir (else executable_dir) 349 return comma-delimited form 350 """ 351 352 if vdir == '': vdir = executable_dir() 353 if vdir == '': return '' 354 355 vpath = '%s/%s' % (vdir, vfile) 356 357 if not os.path.isfile(vpath): return '' 358 359 vdata = read_text_file(vpath, verb=0) 360 if vdata == '': return '' 361 362 return delim.join(vdata) 363 364def write_to_timing_file(data, fname='', nplaces=-1, verb=1): 365 """write the data in stim_times format, over rows 366 (this is not for use with married timing, but for simple times)""" 367 368 if fname == '': return 369 370 fp = open(fname, 'w') 371 if not fp: 372 print("** failed to open '%s' for writing timing" % fname) 373 return 1 374 375 if verb > 0: 376 print("++ writing %d timing rows to %s" % (len(data), fname)) 377 378 fp.write(make_timing_data_string(data, nplaces=nplaces, flag_empty=1, 379 verb=verb)) 380 fp.close() 381 382 return 0 383 384def make_timing_data_string(data, row=-1, nplaces=3, flag_empty=0, 385 mesg='', verb=1): 386 """return a string of row data, to the given number of decimal places 387 if row is non-negative, return a string for the given row, else 388 return a string of all rows""" 389 390 if verb > 2: 391 print('++ make_data_string: row = %d, nplaces = %d, flag_empty = %d' \ 392 % (row, nplaces, flag_empty)) 393 394 if row >= 0: 395 return make_single_row_string(data[row], row, nplaces, flag_empty) 396 397 # make it for all rows 398 if len(mesg) > 0: rstr = "%s :\n" % mesg 399 else: rstr = '' 400 for ind in range(len(data)): 401 rstr += make_single_row_string(data[ind], ind, nplaces, flag_empty) 402 403 return rstr 404 405def make_single_row_string(data, row, nplaces=3, flag_empty=0): 406 """return a string of row data, to the given number of decimal places 407 if row is non-negative, return a string for the given row""" 408 409 rstr = '' 410 411 # if flagging an empty run, use '*' characters 412 if len(data) == 0 and flag_empty: 413 if row == 0: rstr += '* *' 414 else: rstr += '*' 415 416 for val in data: 417 if nplaces >= 0: rstr += '%.*f ' % (nplaces, val) 418 else: rstr += '%g ' % (val) 419 420 return rstr + '\n' 421 422def quotize_list(inlist, opt_prefix='', skip_first=0, quote_wild=0, 423 quote_chars='', ok_chars=''): 424 """given a list of text elements, return a new list where any existing 425 quotes are escaped, and then if there are special characters, put the 426 whole string in single quotes 427 428 if the first character is '-', opt_prefix will be applied 429 if skip_first, do not add initial prefix 430 if quote_wild, quotize any string with '*' or '?', too 431 432 add quote_chars to quote list, remove ok_chars 433 """ 434 if not inlist or len(inlist) < 1: return inlist 435 436 # okay, we haven't yet escaped any existing quotes... 437 438 # default to ignoring wildcards, can always double-nest if needed 439 if quote_wild: qlist = "[({*? " 440 else: qlist = "[({ " 441 442 for c in quote_chars: 443 if not c in qlist: qlist += c 444 for c in ok_chars: 445 posn = qlist.find(c) 446 if posn >= 0: qlist = qlist[0:posn]+qlist[posn+1:] 447 448 newlist = [] 449 first = 1 # ugly, but easier for option processing 450 for qstr in inlist: 451 prefix = '' 452 if skip_first and first: first = 0 # use current (empty) prefix 453 elif len(qstr) == 0: pass 454 elif qstr[0] == '-': prefix = opt_prefix 455 456 quotize = 0 457 for q in qlist: 458 if q in qstr: 459 quotize = 1 460 break 461 if quotize: newlist.append("'%s%s'" % (prefix,qstr)) 462 else: newlist.append("%s%s" % (prefix,qstr)) 463 464 return newlist 465 466def args_as_command(args, prefix='', suffix=''): 467 """given an argument list (such as argv), create a command string, 468 including any prefix and/or suffix strings""" 469 470 if len(args) < 1: return 471 472 cstr = "%s %s" % (os.path.basename(args[0]), 473 ' '.join(quotize_list(args[1:],''))) 474 fstr = add_line_wrappers('%s%s%s' % (prefix,cstr,suffix)) 475 476 return fstr 477 478def show_args_as_command(args, note='command:'): 479 """print the given argument list as a command 480 (this allows users to see wildcard expansions, for example)""" 481 482 print(args_as_command(args, 483 "----------------------------------------------------------------------\n" 484 "%s\n\n " % note, 485 "\n----------------------------------------------------------------------" 486 )) 487 488def exec_tcsh_command(cmd, lines=0, noblank=0, showproc=0): 489 """execute cmd via: tcsh -cf "cmd" 490 return status, output 491 if status == 0, output is stdout 492 else output is stderr+stdout 493 494 if showproc, show command and text output, and do not return any text 495 496 if lines: return a list of lines 497 if noblank (and lines): omit blank lines 498 """ 499 500 # if showproc, show all output immediately 501 if showproc: capture = 0 502 else: capture = 1 503 504 # do not re-process .cshrc, as some actually output text 505 cstr = 'tcsh -cf "%s"' % cmd 506 if showproc: 507 print("%s" % cmd) 508 sys.stdout.flush() 509 status, so, se = BASE.simple_shell_exec(cstr, capture=capture) 510 511 if not status: otext = so 512 else: otext = se+so 513 if lines: 514 slines = otext.splitlines() 515 if noblank: slines = [line for line in slines if line != ''] 516 return status, slines 517 518 return status, otext 519 520def limited_shell_exec(command, nlines=0): 521 """run a simple shell command, returning the top nlines""" 522 st, so, se = BASE.shell_exec2(command, capture=1) 523 if nlines > 0: 524 so = so[0:nlines] 525 se = se[0:nlines] 526 return st, so, se 527 528def write_afni_com_history(fname, length=0, wrap=1): 529 """write the afni_com history to the given file 530 531 if length > 0: limit to that number of entries 532 """ 533 com = BASE.shell_com('hi there') 534 hist = com.shell_history() 535 script = '\n'.join(hist)+'\n' 536 write_text_to_file(fname, script, wrap=wrap) 537 538def get_process_depth(pid=-1, prog=None, fast=1): 539 """print stack of processes up to init""" 540 541 pstack = get_process_stack(pid=pid, fast=fast) 542 543 if prog == None: return len(pstack) 544 545 pids = [pp[0] for pp in pstack if pp[3] == prog] 546 return len(pids) 547 548# get/show_process_stack(), get/show_login_shell() 28 Jun 2013 [rickr] 549def get_process_stack(pid=-1, fast=1, verb=1): 550 """the stack of processes up to init 551 552 return an array of [pid, ppid, user, command] elements 553 --> element 0 should always be init 554 """ 555 def get_pid_index(pids, plist, pid): 556 try: 557 pind = pids.index(pid) 558 except: 559 if verb > 1: 560 print('** GPS pid %s not in pid list:\n' % pid) 561 print('%s' % plist) 562 return -1 563 return pind 564 565 def get_ancestry_indlist(pids, ppids, plist, pid=-1): 566 """return bad status if index() fails""" 567 if pid >= 0: mypid = pid 568 else: mypid = os.getpid() 569 pind = get_pid_index(pids, plist, mypid) 570 if pind < 0: return 1, [] 571 indtree = [pind] 572 while mypid > 1: 573 mypid = ppids[pind] 574 pind = get_pid_index(pids, plist, mypid) 575 if pind < 0: return 1, [] 576 indtree.append(pind) 577 return 0, indtree 578 579 if not fast: 580 stack = get_process_stack_slow(pid=pid) 581 stack.reverse() 582 return stack 583 584 if verb < 2: cmd = 'ps -eo pid,ppid,user,comm' 585 else: cmd = 'ps -eo pid,ppid,user,args' 586 ac = BASE.shell_com(cmd, capture=1) 587 ac.run() 588 if ac.status: 589 print('** GPS command failure for: %s\n' % cmd) 590 print('error output:\n%s' % '\n'.join(ac.se)) 591 return [] 592 593 plist = [ll.split() for ll in ac.so] 594 names = plist[0] 595 plist = plist[1:] 596 597 try: 598 pids = [int(psinfo[0]) for psinfo in plist] 599 ppids = [int(psinfo[1]) for psinfo in plist] 600 except: 601 print('** GPS type failure in plist\n%s' % plist) 602 return [] 603 604 # maybe the ps list is too big of a buffer, so have a backup plan 605 rv, indlist = get_ancestry_indlist(pids, ppids, plist, pid=pid) 606 # if success, set stack, else get it from backup function 607 if rv == 0: stack = [plist[i] for i in indlist] 608 else: stack = get_process_stack_slow(pid=pid) 609 stack.reverse() 610 611 return stack 612 613def get_process_stack_slow(pid=-1, verb=1): 614 """use repeated calls to get stack: 615 ps h -o pid,ppid,user,comm -p PID 616 617 if pid >= 0, use as seed, else use os.getpid() 618 """ 619 if verb > 1: base_cmd = 'ps h -o pid,ppid,user,args -p' 620 else: base_cmd = 'ps h -o pid,ppid,user,comm -p' 621 622 def get_cmd_entries(cmd): 623 ac = BASE.shell_com(cmd, capture=1) 624 ac.run() 625 if ac.status: 626 print('** GPSS command failure for: %s\n' % cmd) 627 print('error output:\n%s' % '\n'.join(ac.se)) 628 return 1, [] 629 ss = ac.so[0] 630 entries = ss.split() 631 if len(entries) == 0: return 1, [] 632 633 return 0, entries 634 635 def get_ppids(cmd, entries): 636 try: 637 pid = int(entries[0]) 638 ppid = int(entries[1]) 639 except: 640 print('** bad GPSS entries for cmd: %s\n %s' % (cmd, entries)) 641 return 1, -1, -1 642 return 0, pid, ppid 643 644 # get pid and ppid 645 if pid >= 0: mypid = pid 646 else: mypid = os.getpid() 647 cmd = '%s %s' % (base_cmd, mypid) 648 rv, entries = get_cmd_entries(cmd) 649 if rv: 650 if mypid == pid: print('** process ID %d not found' % pid) 651 else: print('** getpid() process ID %d not found' % pid) 652 return [] 653 rv, mypid, ppid = get_ppids(cmd, entries) 654 if rv: return [] 655 656 stack = [entries] # entries is valid, so init stack 657 while mypid > 1: 658 cmd = '%s %s' % (base_cmd, ppid) 659 rv, entries = get_cmd_entries(cmd) 660 if rv: return [] 661 rv, mypid, ppid = get_ppids(cmd, entries) 662 if rv: return [] 663 stack.append(entries) 664 665 return stack 666 667def show_process_stack(pid=-1,fast=1,verb=1): 668 """print stack of processes up to init""" 669 pstack = get_process_stack(pid=pid,fast=fast,verb=verb) 670 if len(pstack) == 0: 671 print('** empty process stack') 672 return 673 ulist = [pp[2] for pp in pstack] 674 ml = max_len_in_list(ulist) 675 header = ' PID PPID [USER]' 676 dashes = ' ---- ---- ------' 677 form = '%6s %6s [%s]' 678 ilen = len(form)+4+ml 679 680 print('%-*s : %s' % (ilen, header, 'COMMAND')) 681 print('%-*s %s' % (ilen, dashes, '-------')) 682 for row in pstack: 683 ss = form % (row[0], row[1], row[2]) 684 if len(row) > 3: rv = ' '.join(row[3:]) 685 else: rv = row[3] 686 print('%-*s : %s' % (ilen, ss, rv)) 687 688def get_login_shell(): 689 """return the apparent login shell 690 from get_process_stack(), search down s[3] until a shell is found 691 """ 692 693 pstack = get_process_stack() 694 if len(pstack) == 0: 695 print('** cannot detect shell: empty process stack') 696 return 697 698 # start from init and work down to find first valid shell 699 for pline in pstack: 700 shell = shell_name_strip(pline[3]) 701 if shell: return shell 702 703 return 'SHELL_NOT_DETECTED' 704 705def get_current_shell(): 706 """return the apparent login shell 707 from get_process_stack(), search down s[3] until a shell is found 708 """ 709 710 pstack = get_process_stack() 711 if len(pstack) == 0: 712 print('** cannot detect shell: empty process stack') 713 return 714 pstack.reverse() 715 716 # start from init and work down to find first valid shell 717 for pline in pstack: 718 shell = shell_name_strip(pline[3]) 719 if shell: return shell 720 721 return 'SHELL_NOT_DETECTED' 722 723def shell_name_strip(name): 724 """remove any leading dash or path, then return if valid shell, else ''""" 725 global g_valid_shells 726 if len(name) < 2: return '' 727 728 # strip any leading '-' or path 729 if name[0] == '-': name = name[1:] 730 if name[0] == '/': name = name.split('/')[-1] 731 732 if name in g_valid_shells: return name 733 734 return '' # not a shell 735 736def show_login_shell(verb=0): 737 """print the apparent login shell 738 739 from get_process_stack(), search down s[3] until a shell is found 740 """ 741 shells = ['csh','tcsh','sh','bash','zsh'] 742 743 pstack = get_process_stack() 744 if len(pstack) == 0: 745 print('** cannot detect shell: empty process stack') 746 return 747 748 # start from init and work down to find first valid shell 749 shell = '' 750 for pline in pstack: 751 shell = shell_name_strip(pline[3]) 752 if not shell: continue 753 if verb: print('apparent login shell: %s' % shell) 754 else: print('%s' % shell) 755 break 756 757 if shell == '': 758 if verb: 759 print('** failed to determine login shell, see process stack...\n') 760 show_process_stack() 761 return 762 763 # in verbose mode, see if parent shell is different from login 764 if verb: 765 pstack.reverse() 766 for pline in pstack: 767 sh = shell_name_strip(pline[3]) 768 if not sh: continue 769 770 if sh != shell: print('differs from current shell: %s' % sh) 771 break 772 773def get_unique_sublist(inlist, keep_order=1): 774 """return a copy of inlist, but where elements are unique 775 776 if keep_order, the order is not altered (first one found is kept) 777 (easy to code, but maybe slow) 778 else, sort (if needed), and do a linear pass 779 780 tested with: 781 llist = [3, 4, 7, 4, 5,5,5, 4, 7, 9] 782 get_unique_sublist() 783 get_unique_sublist(keep_order=0) 784 llist.sort() 785 get_unique_sublist(keep_order=0) 786 """ 787 788 if len(inlist) == 0: return [] 789 790 # if keep_order, be slow 791 if keep_order: 792 newlist = [] 793 for val in inlist: 794 if not val in newlist: newlist.append(val) 795 return newlist 796 797 # else, sort only if needed 798 if vals_are_sorted(inlist): 799 slist = inlist 800 else: 801 slist = inlist[:] 802 slist.sort() 803 804 newlist = slist[0:1] 805 for ind in range(len(slist)-1): 806 # look for new values 807 if slist[ind+1] != slist[ind]: newlist.append(slist[ind+1]) 808 809 return newlist 810 811def uniq_list_as_dsets(dsets, byprefix=0, whine=0): 812 """given a list of text elements, create a list of afni_name elements, 813 and check for unique prefixes""" 814 815 if not dsets or len(dsets) < 2: return 1 816 817 if type(dsets[0]) == str: 818 anlist = [BASE.afni_name(dset) for dset in dsets] 819 elif isinstance(dsets[0], BASE.afni_name): 820 anlist = dsets 821 else: 822 print('** ULAD: invalid type for dset list, have value %s' % dsets[0]) 823 return 0 824 825 if byprefix: 826 plist = [an.prefix for an in anlist] 827 else: 828 plist = dsets[:] 829 plist.sort() 830 831 # iterate over dsets, searching for matches 832 uniq = 1 833 for ind in range(len(plist)-1): 834 if anlist[ind].prefix == anlist[ind+1].prefix: 835 uniq = 0 836 break 837 838 if not uniq and whine: 839 print("-----------------------------------------------------------\n" \ 840 "** dataset names are not unique\n\n" \ 841 " (#%d == #%d, '%s' == '%s')\n\n" \ 842 " note: if using a wildcard, please specify a suffix,\n" \ 843 " otherwise datasets may be listed twice\n" \ 844 " e.g. bad use: ED_r*+orig*\n" \ 845 " e.g. good use: ED_r*+orig.HEAD\n" \ 846 "-----------------------------------------------------------\n" \ 847 % (ind, ind+1, anlist[ind].pve(), anlist[ind+1].pve())) 848 849 return uniq 850 851def uniq_list_as_dset_names(dsets, whine=0): 852 """like as_dsets, but the input is a simlpe array of names""" 853 asets = list_to_datasets(dsets, whine=whine) 854 if not asets: return 0 855 return uniq_list_as_dsets(asets, whine=whine) 856 857def list_to_datasets(words, whine=0): 858 """given a list, return the list of afni_name elements 859 - the list can include wildcarding 860 - they must be valid names of existing datasets 861 - return None on error""" 862 863 if not words or len(words) < 1: return [] 864 dsets = [] 865 wlist = [] 866 errs = 0 867 for word in words: 868 glist = glob.glob(word) # first, check for globbing 869 if glist: 870 glist.sort() 871 wlist += glist 872 else: wlist.append(word) 873 # now process all words 874 for word in wlist: 875 dset = BASE.afni_name(word) 876 if dset.exist(): 877 dsets.append(dset) 878 else: 879 if whine: print("** no dataset match for '%s'" % word) 880 errs = 1 881 882 if errs: 883 if whine: print('') # for separation 884 return None 885 return dsets 886 887def dset_prefix_endswith(dname, suffix): 888 """return 1 if an afni_name dset based on dname has a prefix 889 that ends with suffix""" 890 aname = BASE.afni_name(dname) 891 rv = aname.prefix.endswith(suffix) 892 del(aname) 893 return rv 894 895def basis_has_known_response(basis, warn=0): 896 """given a string, if the prefix is either GAM or BLOCK, then the basis 897 function has a known response curve 898 899 if warn, warn users about any basis function peculiarities""" 900 if not basis: return 0 901 902 if starts_with_any_str(basis, basis_known_resp_l): return 1 903 return 0 904 905def basis_is_married(basis): 906 """if the given basis function is known to require married times, return 1 907 """ 908 if not basis: return 0 909 910 if starts_with(basis, 'dmBLOCK'): return 1 911 else: return 0 912 913def basis_has_one_reg(basis, st='times'): 914 """if the given basis function is known to have 1 regressor, return 1 915 """ 916 if not basis: return 0 917 918 # only 'times', 'file' and 'AM1' are acceptable 919 if not st in stim_types_one_reg: return 0 920 921 if starts_with_any_str(basis, basis_one_regr_l): return 1 922 923 return 0 924 925def starts_with(word, sstr): 926 """return 1 if word starts with sstr""" 927 slen = len(sstr) 928 if word[0:slen] == sstr: return 1 929 return 0 930 931def starts_with_any_str(word, slist): 932 """return 1 if word starts with anything in slist""" 933 for sstr in slist: 934 slen = len(sstr) 935 if word[0:slen] == sstr: return 1 936 return 0 937 938def get_default_polort(tr, reps): 939 """compute a default run polort, as done in 3dDeconvolve 940 (leave run_time as TR*reps, rather than TR*(reps-1), to match 3dD) 941 """ 942 943 if tr <= 0 or reps <= 0: 944 print("** cannot guess polort from tr = %f, reps = %d" % (tr,reps)) 945 return 2 # return some default 946 947 return run_time_to_polort(tr*reps) 948 949def run_time_to_polort(run_time): 950 """direct computation: 1+floor(run_time/150)""" 951 return 1+math.floor(run_time/150.0) 952 953def index_to_run_tr(index, rlens, rstyle=1, whine=1): 954 """given index and a list of run lengths, 955 return the corresponding run and TR index 956 957 rstyle: 0/1 whether the run index is 0-based or 1-based 958 959 any negative return indicates an error 960 """ 961 if index < 0: 962 if whine: print('** ind2run_tr: illegal negative index: %d' % index) 963 return 1, index 964 if len(rlens) == 0: 965 if whine: print('** ind2run_tr: missing run lengths') 966 return 1, index 967 968 # if there is only 1 run and it is short, compute modulo 969 rlengths = rlens 970 if len(rlens) == 1 and index >= rlens[0]: 971 rind = index // rlens[0] 972 cind = index % rlens[0] 973 if rstyle: return rind+1, cind 974 else: return rind, cind 975 976 cind = index 977 for rind, nt in enumerate(rlengths): 978 if nt < 0: 979 if whine: 980 print('** ind2run_tr: have negative run length in %s' % rlengths) 981 return -1, -1 982 if cind < nt: 983 if rstyle: return rind+1, cind 984 else: return rind, cind 985 cind -= nt 986 987 if whine: print('** ind2run_tr, index %d outside run list %s' \ 988 % (index, rlengths)) 989 return rind, cind+nt 990 991 992def get_num_warp_pieces(dset, verb=1): 993 """return the number of pieces in the the WARP_DATA transformation 994 for this dataset 995 996 Note that 12 (30 float) pieces would imply manual tlrc, while 1 piece 997 would suggest auto, +acpc, or perhaps an external transform. 998 999 dset : (string) name of afni dataset 1000 1001 if len is a multiple of 30, return len(WARP_DATA)//30 1002 else return 0""" 1003 1004 err, result = get_typed_dset_attr_list(dset, 'WARP_DATA', float, verb=verb) 1005 1006 if err: return 0 # errors printed in called function 1007 1008 nvals = len(result) 1009 npieces = nvals//30 1010 if npieces * 30 != nvals: 1011 print('** GNWP: invalid WARP_DATA length %d' % nvals) 1012 return 0 1013 1014 if verb > 1: print('-- dset %s has a %d-piece warp' % (dset, npieces)) 1015 1016 del(result) 1017 1018 return npieces 1019 1020def get_typed_dset_attr_list(dset, attr, atype, verb=1): 1021 """given an AFNI dataset, return err (0=success), [typed attr list] 1022 1023 dset : (string) name of afni dataset 1024 attr : (string) attribute name 1025 atype : (string) return type of attribute list 1026 verb : verbose level 1027 1028 This ends up running 3dAttribute for the given dataset name. 1029 """ 1030 1031 alist = BASE.read_attribute(dset, attr, verb=verb) 1032 if alist == None and verb > 0: 1033 print("** GTDAL: failed to read dset attr %s, dset = %s" % (attr,dset)) 1034 return 1, [] 1035 1036 err = 0 1037 try: result = [atype(val) for val in alist] 1038 except: 1039 if verb > 0: 1040 print("** GTDAL: failed to convert attr %s to type %s for dset %s"%\ 1041 (attr, atype, dset)) 1042 err = 1 1043 result = [] 1044 1045 return err, result 1046 1047def get_dset_history(dname, lines=1): 1048 """run 3dinfo -history on 'dname' and return the list of lines""" 1049 command = '3dinfo -history %s' % dname 1050 status, olines = exec_tcsh_command(command, lines=lines, noblank=1) 1051 if status: return [] 1052 return olines 1053 1054def get_last_history_command(dname, substr): 1055 """run 3dinfo -history on 'dname' and return the last line 1056 containing 'substr' 1057 return one text line, or '' on failure 1058 """ 1059 olines = get_dset_history(dname) 1060 olen = len(olines) 1061 if olen == 0: return '' 1062 1063 # work backwards 1064 for ind in range(olen-1, -1, -1): 1065 if olines[ind].find(substr) >= 0: return olines[ind] 1066 1067 return '' 1068 1069def get_last_history_ver_pack(dname): 1070 """run 3dinfo -history on 'dname' and return the last line 1071 containing a valid version 1072 return status and the version and package strings 1073 status = 0 on success, 1 on error 1074 """ 1075 olines = get_dset_history(dname) 1076 olen = len(olines) 1077 if olen == 0: return '' 1078 1079 # work backwards and extract the first valid version found 1080 for ind in range(olen-1, -1, -1): 1081 st, vstrs = find_afni_history_version(olines[ind]) 1082 if st == 0: 1083 return 0, vstrs[0], vstrs[1] 1084 1085 return 1, '', '' 1086 1087def get_last_history_version(dname): 1088 """get_last_history_version(DSET_NAME): 1089 return the last AFNI version from the HISTORY 1090 return an empty string on error 1091 """ 1092 st, aver, package = get_last_history_ver_pack(dname) 1093 return aver 1094 1095def get_last_history_package(dname): 1096 """get_last_history_package(DSET_NAME): 1097 return the last AFNI package from the HISTORY 1098 return an empty string on error 1099 """ 1100 st, aver, package = get_last_history_ver_pack(dname) 1101 return package 1102 1103def find_afni_history_version(av_str): 1104 """given {AFNI_A.B.C:PACKAGE}, 1105 return the status and [AFNI_A.B.C, PACKAGE] pair as a list 1106 return 1, [] on error 1107 """ 1108 re_format = '{(AFNI_\d+\.\d+\.\d+):(.*)}' 1109 1110 try: match = re.search(re_format, av_str) 1111 except: return 1, [] 1112 1113 if match is None: 1114 return 1, [] 1115 1116 # we have a match, just verify that we found all 4 pieces 1117 rv = list(match.groups()) 1118 1119 if len(rv) != 2: return 1, [] 1120 1121 return 0, rv 1122 1123def parse_afni_version(av_str): 1124 """given 'AFNI_18.2.10', return status 0 and the int list [18,2,10] 1125 return 1, [] on error 1126 """ 1127 re_format = 'AFNI_(\d+)\.(\d+)\.(\d+)' 1128 1129 try: match = re.search(re_format, av_str) 1130 except: return 1, [] 1131 1132 if match is None: 1133 return 1, [] 1134 1135 # we have a match, convert to ints 1136 rv = list(match.groups()) 1137 1138 if len(rv) != 3: return 1, [] 1139 1140 # return what we have 1141 try: return 0, [int(val) for val in rv] 1142 except: return 1, [] 1143 1144def get_3dinfo(dname, lines=0, verb=0): 1145 """run 3dinfo, possibly -verb or -VERB 1146 if lines, splitlines 1147 return text, or None on failure (applies to string or lines) 1148 """ 1149 vstr = ' ' 1150 if verb == 1: vstr = ' -verb' 1151 elif verb > 1: vstr = ' -VERB' 1152 command = '3dinfo%s %s' % (vstr, dname) 1153 status, output = exec_tcsh_command(command, lines=lines, noblank=1) 1154 if status: return None 1155 1156 return output 1157 1158def get_3dinfo_nt(dname, verb=1): 1159 """run 3dinfo -nt 1160 1161 return 0 on failure (>= 0 on success) 1162 """ 1163 command = '3dinfo -nt %s' % dname 1164 status, output, se = limited_shell_exec(command, nlines=1) 1165 if status or len(output) == 0: 1166 if verb: print('** 3dinfo -nt failure: message is:\n%s%s\n' % (se,output)) 1167 return 0 1168 1169 output = output[0].strip() 1170 if output == 'NO-DSET' : 1171 if verb: print('** 3dinfo -nt: no dataset %s' % dname) 1172 return 0 1173 1174 nt = 0 1175 try: nt = int(output) 1176 except: 1177 if verb: print('** 3dinfo -nt: cannot get NT from %s, for dset %s' \ 1178 % (output, dname)) 1179 return 0 1180 1181 return nt 1182 1183def get_3dinfo_val(dname, val, vtype, verb=1): 1184 """run 3dinfo -val, and convert to vtype (also serves as a test) 1185 1186 return vtype(0) on failure 1187 """ 1188 command = '3dinfo -%s %s' % (val, dname) 1189 status, output, se = limited_shell_exec(command, nlines=1) 1190 if status or len(output) == 0: 1191 if verb: 1192 print('** 3dinfo -%s failure: message is:\n%s%s\n' % (val, se, output)) 1193 return 0 1194 1195 output = output[0].strip() 1196 if output == 'NO-DSET' : 1197 if verb: print('** 3dinfo -%s: no dataset %s' % (val, dname)) 1198 return 0 1199 1200 dval = 0 1201 try: dval = vtype(output) 1202 except: 1203 # allow conversion from float to int as a backup 1204 fail = 0 1205 if vtype == int: 1206 try: 1207 dval = float(output) 1208 dval = vtype(dval) 1209 except: 1210 fail = 1 1211 if verb and fail: 1212 print("** 3dinfo -%s: cannot get val from %s, for dset %s" \ 1213 % (val, output, dname)) 1214 if fail: return vtype(0) 1215 1216 return dval 1217 1218def get_3dinfo_val_list(dname, val, vtype, verb=1): 1219 """run 3dinfo -val, and convert to vtype (also serves as a test) 1220 1221 return None on failure, else a list 1222 """ 1223 command = '3dinfo -%s %s' % (val, dname) 1224 status, output, se = limited_shell_exec(command, nlines=1) 1225 if status or len(output) == 0: 1226 if verb: 1227 print('** 3dinfo -%s failure: message is:\n%s%s\n' % (val, se, output)) 1228 return None 1229 1230 output = output[0].strip() 1231 if output == 'NO-DSET' : 1232 if verb: print('** 3dinfo -%s: no dataset %s' % (val, dname)) 1233 return None 1234 1235 dlist = string_to_type_list(output, vtype) 1236 if dlist == None and verb: 1237 print("** 3dinfo -%s: cannot get val list from %s, for dset %s" \ 1238 % (val, output, dname)) 1239 1240 return dlist 1241 1242def dset_view(dname): 1243 """return the AFNI view for the given dset""" 1244 command = '3dinfo -av_space %s' % dname 1245 status, output = exec_tcsh_command(command) 1246 if status: return '' 1247 return output.replace('\n', '') 1248 1249def get_3d_statpar(dname, vindex, statcode='', verb=0): 1250 """return a single stat param at the given sub-brick index 1251 if statcode, verify 1252 return -1 on error 1253 """ 1254 ilines = get_3dinfo(dname, lines=1, verb=1) 1255 if ilines == None: 1256 print('** failed get_3dinfo(%s)' % dname) 1257 return -1 1258 1259 N = len(ilines) 1260 1261 # find 'At sub-brick #v' line 1262 sstr = 'At sub-brick #%d' % vindex 1263 posn = -1 1264 for ind, line in enumerate(ilines): 1265 posn = line.find(sstr) 1266 if posn >= 0: break 1267 1268 if posn < 0: 1269 print('** 3d_statpar: no %s[%d]' % (dname, vindex)) 1270 return -1 # failure 1271 1272 # check statcode? 1273 lind = ind + 1 1274 if lind >= N: 1275 if verb > 1: print('** 3d_statpar: no space for statpar line') 1276 return -1 1277 1278 sline = ilines[lind] 1279 plist = sline.split() 1280 if statcode: 1281 olist = find_opt_and_params(sline, 'statcode', 2) 1282 if len(olist) < 3: 1283 print('** 3d_statpar: missing expected statcode') 1284 return -1 1285 code = olist[2] 1286 if code[-1] == ';': code = code[:-1] 1287 if code != statcode: 1288 print('** 3d_statpar: statcode %s does not match expected %s'\ 1289 % (code, statcode)) 1290 return -1 1291 if verb > 2: print('-- found %s' % olist) 1292 1293 # now get something like "statpar = 32 x x" 1294 olist = find_opt_and_params(sline, 'statpar', 4) 1295 if len(olist) < 3: 1296 if verb: print('** 3d_statpar: missing expected statpar') 1297 if verb > 2: print(' found %s in %s' % (olist, sline)) 1298 return -1 1299 if verb > 2: print('-- found %s' % olist) 1300 1301 par = -1 1302 try: par = int(olist[2]) 1303 except: 1304 if verb: print('** 3d_statpar: bad stat par[2] in %s' % olist) 1305 return -1 1306 1307 return par 1308 1309def converts_to_type(val, vtype): 1310 """ converts_to_type(val, vtype) : 1311 return whether vtype(val) succeeds""" 1312 rv = 1 1313 try: vtype(val) 1314 except: rv = 0 1315 1316 return rv 1317 1318def find_opt_and_params(text, opt, nopt=0): 1319 """given some text, return the option with that text, as well as 1320 the following 'nopt' parameters (truncated list if not found)""" 1321 tlist = text.split() 1322 1323 if not opt in tlist: return [] 1324 1325 tind = tlist.index(opt) 1326 1327 return tlist[tind:tind+1+nopt] 1328 1329def get_truncated_grid_dim(dset, verb=1): 1330 """return a new (isotropic) grid dimension based on the current grid 1331 - given md = min(DELTAS), return md truncated to 3 significant bits 1332 (first integer this affects is 9->8, then 11->10, etc.) 1333 - return <= 0 on failure 1334 """ 1335 err, dims = get_typed_dset_attr_list(dset, 'DELTA', float) 1336 if err: return -1 1337 if len(dims) < 1: return -1 1338 for ind in range(len(dims)): 1339 dims[ind] = abs(dims[ind]) 1340 md = min(dims) 1341 # changed 2 -> 4 19 Mar 2010 1342 if md >= 4.0: return math.floor(md) 1343 if md <= 0: 1344 print('** failed to get truncated grid dim from %s' % dims) 1345 return 0 1346 1347 return truncate_to_N_bits(md, 3, verb=verb, method='r_then_t') 1348 1349def truncate_to_N_bits(val, bits, verb=1, method='trunc'): 1350 """truncate the real value to most significant N bits 1351 allow for any real val and positive integer bits 1352 1353 method trunc - truncate to 'bits' significant bits 1354 round - round to 'bits' significant bits 1355 r_then_t - round to 2*bits sig bits, then trunc to bits 1356 """ 1357 1358 # allow any real val 1359 if val == 0.0: return 0.0 1360 if val < 0.0: sign, fval = -1, -float(val) 1361 else: sign, fval = 1, float(val) 1362 1363 if verb > 2: print('T2NB: applying sign=%d, fval=%g' % (sign,fval)) 1364 1365 # if r_then_t, start by rounding to 2*bits, then continue to truncate 1366 meth = method 1367 if method == 'r_then_t': 1368 fval = truncate_to_N_bits(val,2*bits,verb,'round') 1369 meth = 'trunc' 1370 1371 if bits <= 0 or type(bits) != type(1): 1372 print("** truncate to N bits: bad bits = ", bits) 1373 return 0.0 1374 1375 # find integer m s.t. 2^(bits-1) <= 2^m * fval < 2^bits 1376 log2 = math.log(2.0) 1377 m = int(math.ceil(bits-1 - math.log(fval)/log2)) 1378 pm = 2**m 1379 1380 # then (round or) truncate to an actual integer in that range 1381 # and divide by 2^m (cannot be r_then_t here) 1382 if meth == 'round': ival = round(pm * fval) 1383 else: ival = math.floor(pm * fval) 1384 retval = sign*float(ival)/pm 1385 1386 if verb > 2: 1387 print('-- T2NB: 2^%d <= 2^%d * %g < 2^%d' % (bits-1,m,fval,bits)) 1388 print(' ival = %g, returning %g' % (ival,retval)) 1389 1390 return retval 1391 1392def test_truncation(top=10.0, bot=0.1, bits=3, e=0.0000001): 1393 """starting at top, repeatedly truncate to bits bits, and subtract e, 1394 while result is greater than bot""" 1395 1396 print('-- truncating from %g down to %g with %d bits' % (top,bot,bits)) 1397 val = top 1398 while val > bot: 1399 trunc = truncate_to_N_bits(val,bits) 1400 print(val, ' -> ', trunc) 1401 val = trunc - e 1402 1403def get_dset_reps_tr(dset, notr=0, verb=1): 1404 """given an AFNI dataset, return err, reps, tr 1405 1406 if notr: use 3dinfo -nt 1407 1408 err = error code (0 = success, else failure) 1409 reps = number of TRs in the dataset 1410 tr = length of TR, in seconds 1411 """ 1412 1413 # use 3dinfo directly, instead of TAXIS attributes 30 Jun 2016 1414 1415 reps = get_3dinfo_val(dset, 'nt', int, verb=verb) 1416 tr = get_3dinfo_val(dset, 'tr', float, verb=verb) 1417 1418 # store timing info in a list (to get reps and timing units) 1419 if reps == 0: 1420 print("** failed to find the number of TRs from dset '%s'" % dset) 1421 return 1, None, None 1422 1423 if verb > 1: print('-- dset %s : reps = %d, tr = %ss' % (dset, reps, tr)) 1424 1425 return 0, reps, tr 1426 1427def gaussian_width_to_fwhm(width, mode): 1428 """convert the given 'width' of gaussian 'mode' to FWHM 1429 (full width at half max) 1430 1431 mode can be one of: fwhm, rms, sigma 1432 1433 conversion based on valid mode is: 1434 rms = 0.42466090 * fwhm 1435 sigma = 0.57735027 * rms 1436 1437 implying: 1438 fwhm = 2.354820 * rms 1439 fwhm = 4.078668 * sigma 1440 1441 return 0 on failure or error""" 1442 1443 if width <= 0.0: return 0.0 1444 if mode == 'fwhm': return width 1445 if mode == 'rms': return width * 2.354820 1446 if mode == 'sigma': return width * 4.078668 1447 1448 print("** GW2F: illegal mode '%s'" % mode) 1449 1450 return 0.0 1451 1452def attr_equals_val(object, attr, val): 1453 """return 1 of the object has attribute attr and it equals val""" 1454 1455 rv = 0 1456 try: 1457 oval = getattr(object, attr) 1458 if oval == val: rv = 1 1459 except: pass 1460 1461 return rv 1462 1463# ---------------------------------------------------------------------- 1464# begin matrix functions 1465 1466def num_cols_1D(filename): 1467 """return the number of columns in a 1D file""" 1468 mat = TD.read_1D_file(filename) 1469 if not mat or len(mat) == 0: return 0 1470 return len(mat[0]) 1471 1472def num_rows_1D(filename): 1473 """return the number of columns in a 1D file""" 1474 mat = TD.read_1D_file(filename) 1475 if not mat: return 0 1476 return len(mat) 1477 1478def max_dim_1D(filename): 1479 """return the larger of the number of rows or columns""" 1480 mat = TD.read_1D_file(filename) 1481 if not mat: return 0 1482 rows = len(mat) 1483 cols = len(mat[0]) 1484 if rows >= cols: return rows 1485 else: return cols 1486 1487def is_matrix_square( mat, full_check=False ): 1488 """return {0|1} about matrix being square""" 1489 1490 if not(mat) or len(mat)==0 : return 0 1491 1492 rows = len(mat) 1493 1494 # can do complete check, or just [0]th ele check 1495 N = 1 1496 if full_check : 1497 N = rows 1498 1499 for i in range(N): 1500 if len(mat[i]) != rows : 1501 return 0 1502 1503 # have we survived to here? then -> square. 1504 return 1 1505 1506def mat_row_mincol_maxcol_ragged_square(M): 1507 """check 5 properties of a matrix (list of lists) and return 5 ints: 1508 nrow 1509 min ncol 1510 max ncol 1511 is_ragged 1512 is_square 1513 1514 """ 1515 1516 if not(M): return 0,0,0,0,0 1517 1518 is_square = 0 # just default; can change below 1519 1520 nrow = len(M) 1521 all_clen = [len(r) for r in M] 1522 ncolmin = min(all_clen) 1523 ncolmax = max(all_clen) 1524 1525 if ncolmin == ncolmax : 1526 is_ragged = 0 1527 if ncolmin == nrow : 1528 is_square = 1 1529 else: 1530 is_ragged = 1 1531 1532 return nrow, ncolmin, ncolmax, is_ragged, is_square 1533 1534def transpose(matrix): 1535 """transpose a 2D matrix, returning the new one""" 1536 if matrix is None: return [] 1537 1538 rows = len(matrix) 1539 1540 # handle trivial case 1541 if rows == 0: return [] 1542 1543 cols = len(matrix[0]) 1544 newmat = [] 1545 for c in range(cols): 1546 newrow = [] 1547 for r in range(rows): 1548 newrow.append(matrix[r][c]) 1549 newmat.append(newrow) 1550 return newmat 1551 1552def derivative(vector, in_place=0, direct=0): 1553 """take the derivative of the vector, setting d(t=0) = 0 1554 1555 in_place: if set, the passed vector will be modified 1556 direct: if 0, use backward difference, if 1 use forward 1557 1558 return v[t]-v[t-1] for 1 in {1,..,len(v)-1}, d[0]=0""" 1559 1560 if type(vector) != type([]): 1561 print("** cannot take derivative of non-vector '%s'" % vector) 1562 return None 1563 1564 if in_place: vec = vector # reference original 1565 else: vec = vector[:] # start with copy 1566 1567 # count from the end to allow memory overwrite 1568 if direct: # forward difference 1569 vlen = len(vec) 1570 for t in range(vlen-1): 1571 vec[t] = vec[t+1] - vec[t] 1572 vec[vlen-1] = 0 1573 else: # backward difference 1574 for t in range(len(vec)-1, 0, -1): 1575 vec[t] -= vec[t-1] 1576 vec[0] = 0 1577 1578 return vec 1579 1580def matrix_multiply_2D(A, B, zero_dtype=''): 1581 '''Perform matrix multiplication of 2D lists, A and B, which can be of 1582arbitrary dimension (subject to Arow/Bcol matching). 1583 1584 Each input must have 2 indices. So, the following are valid inputs: 1585 x = [[1, 2, 3]] 1586 while these are NOT valid: 1587 x = [1, 2, 3] 1588 1589 Output a new array of appropriate dims, which will also always 1590 have 2 inds (unless one input had zero row or col; in this case, 1591 return empty arr). So, output might look like: 1592 [[16], [10], [12]] 1593 [[2]] 1594 [] 1595 etc. 1596 1597 The dtype of zeros used to initialize the matrix will be either: 1598 + whatever the [0][0]th element of the A matrix is, or 1599 + whatever is specified with 'zero_dtype' (which can be: bool, 1600 int, float, complex, ...) 1601 NB: any valid zero_dtype value will take precedence. 1602 1603 ''' 1604 1605 # ------------ check dims/shapes ------------ 1606 NrowA = len(A) 1607 if not(NrowA) : 1608 print("** ERROR: No rows in matrix A\n") 1609 sys.exit(4) 1610 NcolA = len(A[0]) 1611 NrowB = len(B) 1612 if not(NrowB) : 1613 print("** ERROR: No rows in matrix B\n") 1614 sys.exit(4) 1615 NcolB = len(B[0]) 1616 1617 # Zero row/col case: return empty 1D arr 1618 if not(NcolA) or not(NcolB) or not(NrowA) or not(NrowB): 1619 return [] 1620 1621 if NcolA != NrowB : 1622 print("** ERROR: Ncol in matrix A ({}) != Nrow in matrix B ({})." 1623 "".format(NcolA, NrowB)) 1624 sys.exit(5) 1625 1626 # Allow for different types of matrices 1627 ZZ = calc_zero_dtype(A[0][0], zero_dtype) 1628 1629 # Initialize output list of correct dimensions 1630 C = [[ZZ] * NcolB for row in range(NrowA)] 1631 1632 # Ye olde matrix multiplication, boolean style 1633 if type(ZZ) == bool : 1634 for i in range(NrowA): 1635 for j in range(NcolB): 1636 for k in range(NcolA): 1637 C[i][j] = C[i][j] or (A[i][k] and B[k][j]) 1638 else: 1639 # Ye olde matrix multiplication 1640 for i in range(NrowA): 1641 for j in range(NcolB): 1642 for k in range(NcolA): 1643 C[i][j]+= A[i][k] * B[k][j] 1644 1645 return C 1646 1647def matrix_sum_abs_val_ele_row(A, zero_dtype='' ): 1648 '''Input is a 2D list (matrix). 1649 1650 Output is a list of sums of the absolute values of elements in 1651 rows (SAVERs!). 1652 1653 ''' 1654 1655 N = len(A) 1656 if not(N) : 1657 return [] 1658 Ncol = len(A[0]) # if this doesn't exist, someone will be unhappy 1659 1660 # initialize for summation 1661 ZZ = calc_zero_dtype(A[0][0], zero_dtype) 1662 out = [ZZ]*N 1663 1664 # sum abs vals across rows 1665 if type(ZZ) == bool : 1666 for i in range(N): 1667 for j in range(len(A[i])) : 1668 out[i] = out[i] or abs(A[i][j]) 1669 else: 1670 for i in range(N): 1671 for j in range(len(A[i])) : 1672 out[i]+= abs(A[i][j]) 1673 1674 return out 1675 1676def calc_zero_dtype(x, zero_dtype=''): 1677 '''Calc an appropriate zero for initializing a list-matrix. 1678 1679 ''' 1680 1681 # priority to an entered zero_dtype 1682 if zero_dtype == int : return 0 1683 elif zero_dtype == float : return 0.0 1684 elif zero_dtype == bool : return False 1685 elif zero_dtype == complex : return 0j 1686 elif type(x) == int : return 0 1687 elif type(x) == float : return 0.0 1688 elif type(x) == bool : return False 1689 elif type(x) == complex : return 0j 1690 else: 1691 print("+* Warning: unrecognized numeric type: {}" 1692 "".format(type(x))) 1693 print(" Initial matrix zeros will be type 'int'") 1694 return 0 1695 1696def make_timing_string(data, nruns, tr, invert=0): 1697 """evaluating the data array as boolean (zero or non-zero), return 1698 non-zero entries in stim_times format 1699 1700 data : single vector (length must be multiple of nruns) 1701 nruns : number of runs (must divide len(data)) 1702 tr : TR in seconds, with data viewed as per TR 1703 1704 invert : (optional) invert the boolean logic (write 0 times) 1705 1706 return err code (0 on success), stim_times string""" 1707 1708 if not data: 1709 print("** make_timing_string: missing data") 1710 return 1, '' 1711 if not type(data) == type([]): 1712 print("** make_timing_string: data is not a list") 1713 return 1, '' 1714 1715 nvals = len(data) 1716 rlen = nvals // nruns 1717 1718 if nruns * rlen != nvals: 1719 print("** make_timing_str: nruns %d does not divide nvals %d"%(rlen,nvals)) 1720 return 1, '' 1721 if tr <= 0.0: 1722 print("** make_timing_string: bad tr = %g" % tr) 1723 return 1, '' 1724 1725 rstr = '' 1726 1727 for run in range(nruns): 1728 bot = run*rlen 1729 if invert: rvals = [1*(data[i] == 0) for i in range(bot,bot+rlen)] 1730 else: rvals = [1*(data[i] != 0) for i in range(bot,bot+rlen)] 1731 # if the run is empty, print 1 or 2 '*' 1732 nzero = rvals.count(0) 1733 if nzero == rlen: 1734 if run == 0: rstr += '* *' 1735 else: rstr += '*' 1736 else: 1737 rstr += ' '.join(['%g'%(i*tr) for i in range(rlen) if rvals[i]]) 1738 1739 # if run0 and exactly 1 non-zero value, print a trailing '*' 1740 if run == 0 and nzero == rlen-1: rstr += ' *' 1741 rstr += '\n' 1742 1743 return 0, rstr 1744 1745def make_CENSORTR_string(data, nruns=0, rlens=[], invert=0, asopt=0, verb=1): 1746 """evaluating the data array as boolean (zero or non-zero), return 1747 non-zero entries in CENSORTR format 1748 1749 data : single vector (length must be multiple of nruns) 1750 nruns : number of runs (must divide len(data)) 1751 (ignored if rlens is provided) 1752 rlens : run lengths (required if run lengths vary) 1753 asopt : if set, return as complete -CENSORTR option 1754 1755 invert : (optional) invert the boolean logic (write 0 TRs) 1756 1757 return err code (0 on success), CENSORTR string""" 1758 1759 if not data: 1760 print("** CENSORTR_str: missing data") 1761 return 1, '' 1762 if not type(data) == type([]): 1763 print("** CENSORTR_str: data is not a list") 1764 return 1, '' 1765 1766 nvals = len(data) 1767 1768 # we must have either nruns or a valid rlens list 1769 if nruns <= 0 and len(rlens) < 1: 1770 print('** make_CENSORTR_string: neither nruns nor rlens') 1771 return 1, '' 1772 1773 if rlens: 1774 rlist = rlens 1775 runs = len(rlist) 1776 else: 1777 rlist = [(nvals//nruns) for run in range(nruns)] 1778 runs = nruns 1779 1780 if verb > 1: 1781 print('-- CENSORTR: applying run lengths (%d) : %s' % (runs, rlist)) 1782 1783 if loc_sum(rlist) != nvals: 1784 print("** CENSORTR_str: sum of run lengths %d != nvals %d" \ 1785 % (loc_sum(rlist),nvals)) 1786 return 1, '' 1787 1788 rstr = '' 1789 1790 bot = 0 1791 for run in range(runs): 1792 rlen = rlist[run] 1793 if invert: rvals = [1*(data[i] == 0) for i in range(bot,bot+rlen)] 1794 else: rvals = [1*(data[i] != 0) for i in range(bot,bot+rlen)] 1795 bot += rlen # adjust bottom index for next run 1796 1797 # if the run is empty, print 1 or 2 '*' 1798 nzero = rvals.count(0) 1799 if nzero == rlen: continue 1800 1801 # make a ',' and '..' string listing TR indices 1802 estr = encode_1D_ints([i for i in range(rlen) if rvals[i]]) 1803 1804 # every ',' separated piece needs to be preceeded by RUN: 1805 rstr += "%d:%s " % (run+1, estr.replace(',', ',%d:'%(run+1))) 1806 1807 if asopt and rstr != '': rstr = "-CENSORTR %s" % rstr 1808 1809 return 0, rstr 1810 1811# ------------------ 1812# funcs for specific kind of matrix: list of 2D list matrices 1813 1814def check_list_2dmat_and_mask(L, mask=None): 1815 """A check to see if L is a list of 2D matrices. 1816 1817 Also, check if any input mask matches the dimensions of the 1818 submatrices of L. 1819 1820 This does not check if all submatrices have the same dimensions. 1821 Maybe someday it will. 1822 1823 Returns 4 integers: 1824 + 'success' value :1 if L is 2D mat and any mask matches 1825 + len(L) 1826 + nrow in L submatrix 1827 + ncol in L submatrix 1828 1829 """ 1830 1831 # calc recursive list of embedded list lengths 1832 Ldims = get_list_mat_dims(L) 1833 1834 # need a [N, nrow, ncol] here 1835 if len(Ldims) != 3 : 1836 BASE.EP("Matrix fails test for being a list of 2D matrices;\m" 1837 "instead of having 3 dims, it has {}".format(len(Ldims))) 1838 1839 if mask != None : 1840 mdims = get_list_mat_dims(mask) 1841 if mdims[0] != Ldims[1] or mdims[1] != Ldims[2] : 1842 BASE.EP("matrix dimensions don't match:\n" 1843 "mask dims: {}, {}" 1844 "each matrix in list has dims: {}, {}" 1845 "".format(mdims[0], mdims[1], Ldims[1], Ldims[2])) 1846 1847 return 1, Ldims[0], Ldims[1], Ldims[2] 1848 1849def calc_list_2dmat_count_nonzero(L, mask=None, mode='count'): 1850 """L is a list of 2D list matrices. Each submatrix must be the same 1851 shape. 1852 1853 Calculate the number of nonzero entries across len(L) for each 1854 submatrix element. 1855 1856 Return a matrix of a different style, depending on the specified 1857 'mode': 1858 'count' :integer number of nonzero elements, i.e., max value is 1859 len(L); this is the default 1860 'frac' :output is floating point: number of nonzero values at a 1861 given matrix element, divided by len(L); therefore, 1862 values here are in a range [0, 1] 1863 'all_nz' :output is a matrix of binary vals, where 1 means that 1864 all values across the length of L are nonzero for that 1865 element, and 0 means that one or more value is zero 1866 'any_nz' :output is a matrix of binary vals, where 1 means that 1867 at least one value across the length of L is nonzero for 1868 that element, and 0 means that all values are zero 1869 1870 A mask can be entered, so that all the above calcs are only done 1871 in the mask; outside of it, all values are 0 or False. 1872 1873 if len(L) == 0: 1874 1 empty matrix is returned 1875 else: 1876 normal calcs 1877 1878 """ 1879 1880 # check if L is a list of 2D mats, and if the mask matches 1881 # submatrices (if used) 1882 IS_OK, N, nrow, ncol = check_list_2dmat_and_mask(L, mask) 1883 1884 mmat = [[0]*ncol for i in range(nrow)] # [PT: June 9, 2020] fixed 1885 Nfl = float(N) 1886 1887 if mask == None : 1888 mask = [[1]*ncol for i in range(nrow)] # [PT: June 9, 2020] fixed 1889 1890 for ii in range(nrow) : 1891 for jj in range(ncol) : 1892 if mask[ii][jj] : 1893 for ll in range(N) : 1894 if L[ll][ii][jj] : 1895 mmat[ii][jj]+= 1 1896 1897 if mode == 'frac': 1898 mmat[ii][jj]/= Nfl 1899 elif mode == 'any_nz' : 1900 mmat[ii][jj] = int(bool(mmat[ii][jj])) 1901 elif mode == 'all_nz' : 1902 if mmat[ii][jj] == N : 1903 mmat[ii][jj] = 1 1904 else: 1905 mmat[ii][jj] = 0 1906 1907 return mmat 1908 1909def calc_list_2dmat_mean_stdev_max_min(L, mask=None, ddof=1 ): 1910 """L is a list of 2D list matrices. Each submatrix must be the same 1911 shape. 1912 1913 Calculate four elementwise properties across the length of L, and 1914 return four separate matrices. 1915 1916 NB: by "list of 2D matrices," we mean that having L[0] = [[1,2,3]] 1917 would lead to actual calcs, but having L[0] = [1,2,3] would not. 1918 That is, each element of L must have both a nonzero number of rows 1919 and cols. 1920 1921 For any shape of list matrix that isn't a "list of 2D matrices", 1922 we return all empty matrices. 1923 1924 mask :An optional mask can be be input; it will be treated as 1925 boolean. Values outside the mask will simply be 0. 1926 1927 ddof :control the 'delta degrees of freedom' in the divisor of 1928 stdev calc; numpy's np.std() has ddof=1 by default, so the 1929 denominator is "N"; here, default is ddof=1, so that 1930 the denominator is N-1 1931 1932 if len(L) == 0: 1933 4 empty matrices are returned 1934 elif len(L) == 1: 1935 mean, min and max matrices are copies of input, and stdev is all 0s 1936 else: 1937 normal calcs 1938 1939 """ 1940 1941 # check if L is a list of 2D mats, and if the mask matches 1942 # submatrices (if used) 1943 IS_OK, N, nrow, ncol = check_list_2dmat_and_mask(L, mask) 1944 1945 if N == 1 : 1946 mmean = copy.deepcopy(L[0]) 1947 mstdev = [[0.]*ncol for i in range(nrow)] # [PT: June 9, 2020] fixed 1948 mmax = copy.deepcopy(L[0]) 1949 mmin = copy.deepcopy(L[0]) 1950 if mask != None : 1951 for ii in range(nrow) : 1952 for jj in range(ncol) : 1953 if not(mask[ii][jj]) : 1954 mmean[ii][jj] = 0 1955 mmin[ii][jj] = 0 1956 mmax[ii][jj] = 0 1957 return mmean, mstdev, mmax, mmin 1958 1959 Nfl = float(N) 1960 Nstdev = Nfl - ddof # ... and we know Nfl>1.0, if here 1961 1962 # Initialize mats 1963 mmean = [[0]*ncol for i in range(nrow)] # [PT: June 9, 2020] fixed 1964 mstdev = [[0]*ncol for i in range(nrow)] # [PT: June 9, 2020] fixed 1965 mmax = copy.deepcopy(L[0]) 1966 mmin = copy.deepcopy(L[0]) 1967 1968 if mask == None : 1969 mask = [[1]*ncol for i in range(nrow)] # [PT: June 9, 2020] fixed 1970 1971 for ii in range(nrow) : 1972 for jj in range(ncol) : 1973 if mask[ii][jj] : 1974 for ll in range(N) : 1975 x = L[ll][ii][jj] 1976 mmean[ii][jj]+= x 1977 mstdev[ii][jj]+= x*x 1978 if x < mmin[ii][jj] : 1979 mmin[ii][jj] = x 1980 if x > mmax[ii][jj] : 1981 mmax[ii][jj] = x 1982 print(mmean[ii][jj],Nfl) 1983 mmean[ii][jj]/= Nfl 1984 mstdev[ii][jj]-= N*mmean[ii][jj]*mmean[ii][jj] 1985 mstdev[ii][jj] = (mstdev[ii][jj]/Nstdev)**0.5 1986 1987 return mmean, mstdev, mmax, mmin 1988 1989def get_list_mat_dims(L, lengths=[]) : 1990 """Calc dims of list-style matrix. Recursive. 1991 1992 Returns a list of lengths of embedded lists; does not check that 1993 each list is the same length at a given level. 1994 1995 """ 1996 1997 if type(L) == list: 1998 a = copy.deepcopy(lengths) 1999 a.append(len(L)) 2000 return get_list_mat_dims(L[0], lengths=a) 2001 else: 2002 return lengths 2003 2004 2005# end matrix functions 2006# ---------------------------------------------------------------------- 2007# index encode/decode functions 2008 2009def encode_1D_ints(ilist): 2010 """convert a list of integers to a ',' and '..' separated string""" 2011 if not ilist: return '' 2012 if len(ilist) < 1: return '' 2013 2014 text = '%d' % ilist[0] 2015 prev = ilist[0] 2016 ind = 1 2017 while ind < len(ilist): 2018 ncontinue = consec_len(ilist, ind-1) - 1 2019 if ncontinue <= 1: # then no '..' continuation, use ',' 2020 text = text + ',%d' % ilist[ind] 2021 ind += 1 2022 else: 2023 text = text + '..%d' % ilist[ind+ncontinue-1] 2024 ind += ncontinue 2025 2026 return text 2027 2028def consec_len(ilist, start): 2029 """return the length of consecutive integers - always at least 1""" 2030 prev = ilist[start] 2031 length = len(ilist) 2032 for ind in range(start+1, length+1): 2033 if ind == length: break 2034 if ilist[ind] != prev + 1: 2035 break 2036 prev = ilist[ind] 2037 if ind == start: length = 1 2038 else: length = ind-start 2039 2040 return length 2041 2042def restrict_by_index_lists(dlist, ilist, base=0, nonempty=1, verb=1): 2043 """restrict elements of dlist by indices in ilist 2044 2045 ilist : can be string or list of strings 2046 (require unique composite list) 2047 base : can be 0 or 1 (0-based or 1-based) 2048 nonempty : if set, sub-lists are not allowed to be empty 2049 verb : verbose level, default is to only report errors 2050 2051 return status, sub-list 2052 status = 0 on success, 1 on error 2053 """ 2054 2055 # if either object is empty, there is nothing to do 2056 if not ilist or not dlist: return 0, [] 2057 2058 if type(ilist) == str: ilist = [ilist] 2059 2060 if base not in [0,1]: 2061 if verb: print('** restrict_by_index_list: bad base = %d' % base) 2062 return 1, [] 2063 2064 # set imax to correctly imply '$' index 2065 if base: imax = len(dlist) # 1-based 2066 else: imax = len(dlist)-1 # 0-based 2067 2068 composite = [] 2069 for ind, istr in enumerate(ilist): 2070 if type(istr) != str: 2071 print('** RBIL: bad index selector %s' % istr) 2072 return 1, [] 2073 curlist = decode_1D_ints(istr, verb=verb, imax=imax) 2074 if not curlist and nonempty: 2075 if verb: print("** empty index list for istr[%d]='%s'" % (ind,istr)) 2076 return 1, [] 2077 composite.extend(curlist) 2078 if verb > 3: print('-- index %d, ilist %s' % (ind, curlist)) 2079 2080 if not vals_are_unique(composite): 2081 if verb: print('** RBIL: composite index list elements are not unique') 2082 return 1, [] 2083 2084 cmin = min(composite) 2085 cmax = max(composite) 2086 if cmin < 0: 2087 if verb: print('** RBIL: cannot choose negative indices') 2088 return 1, [] 2089 elif base and cmin == 0: 2090 if verb: print('** RBIL: 1-based index list seems 0-based') 2091 return 1, [] 2092 elif cmax > imax: 2093 if verb: print('** RBIL: index value %d exceeds %d-based limit %d' \ 2094 % (cmax, base, imax)) 2095 return 1, [] 2096 2097 # now convert 1-based to 0-based, if needed 2098 if base: clist = [v-1 for v in composite] 2099 else: clist = composite 2100 2101 # the big finish 2102 return 0, [dlist[ind] for ind in clist] 2103 2104def decode_1D_ints(istr, verb=1, imax=-1, labels=[]): 2105 """Decode a comma-delimited string of ints, ranges and A@B syntax, 2106 and AFNI-style sub-brick selectors (including A..B(C)). 2107 If the A..B format is used, and B=='$', then B gets 'imax'. 2108 If the list is enclosed in [], <> or ##, strip those characters. 2109 - return a list of ints""" 2110 2111 newstr = strip_list_brackets(istr, verb) 2112 slist = newstr.split(',') 2113 if len(slist) == 0: 2114 if verb > 1: print("-- empty 1D_ints from string '%s'" % istr) 2115 return [] 2116 elif verb > 3: 2117 print("-- decoding stripped list '%s' (nlabs = %d)" \ 2118 % (newstr,len(labels))) 2119 ilist = [] # init return list 2120 for s in slist: 2121 try: 2122 if s.find('@') >= 0: # then expect "A@B" 2123 [N, val] = [n for n in s.split('@')] 2124 N = int(N) 2125 val = to_int_special(val, '$', imax, labels) 2126 ilist.extend([val for i in range(N)]) 2127 elif s.find('..') >= 0: # then expect "A..B" 2128 pos = s.find('..') 2129 if s.find('(', pos) > 0: # look for "A..B(C)" 2130 [v1, v2] = [n for n in s.split('..')] 2131 v1 = to_int_special(v1, '$', imax, labels) 2132 [v2, step] = v2.split('(') 2133 v2 = to_int_special(v2, '$', imax, labels) 2134 # have start and end values, get step 2135 step, junk = step.split(')') 2136 step = int(step) 2137 if step > 0: inc = 1 2138 elif step < 0: inc = -1 2139 else: 2140 print("** decode: illegal step of 0 in '%s'" % istr) 2141 return [] 2142 ilist.extend([i for i in range(v1, v2+inc, step)]) 2143 else: 2144 [v1, v2] = [n for n in s.split('..')] 2145 v1 = to_int_special(v1, '$', imax, labels) 2146 v2 = to_int_special(v2, '$', imax, labels) 2147 if v1 < v2 : step = 1 2148 else: step = -1 2149 ilist.extend([i for i in range(v1, v2+step, step)]) 2150 else: 2151 ilist.extend([to_int_special(s, '$', imax, labels)]) 2152 except: 2153 print("** cannot decode_1D '%s' in '%s'" % (s, istr)) 2154 return [] 2155 if verb > 3: print('++ ilist: %s' % ilist) 2156 del(newstr) 2157 return ilist 2158 2159def to_int_special(cval, spec, sint, labels=[]): 2160 """basically return int(cval), but if cval==spec, return sint 2161 2162 if labels were passed, allow cval to be one 2163 2164 cval: int as character string, or a label 2165 spec: special value as string 2166 sint: special value as int""" 2167 if cval == spec: 2168 return sint 2169 elif cval[0].isalpha() and cval in labels: 2170 return labels.index(cval) 2171 else: 2172 return int(cval) 2173 2174def extract_subbrick_selection(sstring): 2175 """search sstring for something like: [DIGITS*($|DIGITS)] 2176 - do not consider all DIGITS, '..', ',', '(DIGITS)' pieces, 2177 just let '*' refer to anything but another '[' 2178 """ 2179 import re 2180 res = re.search('\[\d+[^\[]*]', sstring) 2181 if res == None: return '' 2182 return res.group(0) 2183 2184def strip_list_brackets(istr, verb=1): 2185 """strip of any [], {}, <> or ## surrounding this string 2186 - can only be surrounded by whitespace 2187 - assume only one pair 2188 - allow the trailing character to be missing 2189 return the remaining string""" 2190 2191 if len(istr) < 1: return istr 2192 2193 icopy = istr.strip() 2194 2195 # strip any of these pairs 2196 for pairs in [ ['[',']'], ['{','}'], ['<','>'], ['#','#'] ]: 2197 if icopy[0] == pairs[0]: 2198 ind1 = icopy.find(pairs[1], 1) 2199 if verb > 1: print('-- stripping %s%s at %d,%d in %s' % \ 2200 (pairs[0],pairs[1],0,ind1,icopy)) 2201 if ind1 > 0: return icopy[1:ind1] 2202 else: return icopy[1:] 2203 2204 if verb > 2: print("-- nothing to strip from '%s'" % icopy) 2205 2206 return icopy 2207 2208def replace_n_squeeze(instr, oldstr, newstr): 2209 """like string.replace(), but remove all spaces around oldstr 2210 (so probably want some space in newstr)""" 2211 # while oldstr is found 2212 # find last preceeding keep posn (before oldstr and spaces) 2213 # find next following keep posn (after oldstr and spaces) 2214 # set result = result[0:first] + newstr + result[last:] 2215 newlen = len(newstr) 2216 result = instr 2217 posn = result.find(oldstr) 2218 while posn >= 0: 2219 rlen = len(result) 2220 start = posn-1 2221 while start >= 0 and result[start] == ' ': start -= 1 2222 if start >= 0: newres = result[0:start+1] + newstr 2223 else: newres = newstr 2224 end = posn + newlen 2225 while end < rlen and result[end] == ' ': end += 1 2226 if end < rlen: newres += result[end:] 2227 2228 result = newres 2229 posn = result.find(oldstr) 2230 2231 return result 2232 2233# ---------------------------------------------------------------------- 2234# line wrapper functions 2235 2236def list_to_wrapped_command(cname, llist, nindent=10, nextra=3, maxlen=76): 2237 """return a string that is a 'cname' command, indented by 2238 nindent, with nextra indentation for option continuation 2239 2240 This function taks a command and a list of options with parameters, 2241 and furnishes a wrapped command, where each option entry is on its 2242 own line, and any option entry line wraps includes nextra indentation. 2243 2244 - wrap each option line without indentation 2245 - split across \\\n; ==> have good indentation 2246 - keep all in array 2247 2248 - join array with indentation and \\\n 2249 - finally, align the \\\n wrappers 2250 """ 2251 ntot = nindent + nextra 2252 isind = ' '*nindent 2253 istot = ' '*ntot 2254 jstr = ' \\\n' + istot 2255 2256 clist = [cname] 2257 for line in llist: 2258 wline = add_line_wrappers(line, maxlen=(maxlen-ntot)) 2259 wlist = wline.split('\\\n') 2260 clist.extend(wlist) 2261 2262 return align_wrappers(isind + jstr.join(clist)) 2263 2264 2265# MAIN wrapper: add line wrappers ('\'), and align them all 2266def add_line_wrappers(commands, wrapstr='\\\n', maxlen=78, verb=1): 2267 """wrap long lines with 'wrapstr' (probably '\\\n' or just '\n') 2268 if '\\\n', align all wrapstr strings""" 2269 new_cmd = '' 2270 posn = 0 2271 2272 while needs_wrapper(commands, maxlen, posn): 2273 2274 end = find_command_end(commands, posn) 2275 2276 if not needs_wrapper(commands, maxlen, posn, end): # command is okay 2277 if end < 0: new_cmd = new_cmd + commands[posn:] 2278 else : new_cmd = new_cmd + commands[posn:end+1] 2279 posn = end+1 2280 continue 2281 2282 # command needs wrapping 2283 new_cmd += insert_wrappers(commands, posn, end, wstring=wrapstr, 2284 maxlen=maxlen, verb=verb) 2285 2286 posn = end + 1 # else, update posn and continue 2287 2288 result = new_cmd + commands[posn:] 2289 2290 # wrappers are in, now align them 2291 if wrapstr == '\\\n': return align_wrappers(result) 2292 else: return result 2293 2294def align_wrappers(command): 2295 """align all '\\\n' strings to be the largest offset 2296 from the previous '\n'""" 2297 2298 # first, find the maximum offset 2299 posn = 0 2300 max = -1 2301 while 1: 2302 next = command.find('\n',posn) 2303 if next < 0: break 2304 if next > posn and command[next-1] == '\\': # check against max 2305 width = next - 1 - posn 2306 if width > max: max = width 2307 posn = next + 1 # look past it 2308 2309 if max < 0: return command # none found 2310 2311 # repeat the previous loop, but adding appropriate spaces 2312 new_cmd = '' 2313 posn = 0 2314 while 1: 2315 next = command.find('\n',posn) 2316 if next < 0: break 2317 if next > posn and command[next-1] == '\\': # check against max 2318 width = next - 1 - posn 2319 if width < max: # then insert appropriate spaces 2320 new_cmd += command[posn:next-1] + ' '*(max-width) + '\\\n' 2321 posn = next + 1 2322 continue 2323 2324 # just duplicate from the previous posn 2325 new_cmd += command[posn:next+1] 2326 posn = next + 1 # look past it 2327 2328 if posn < len(command): new_cmd += command[posn:] 2329 2330 return new_cmd 2331 2332def insert_wrappers(command, start=0, end=-1, wstring='\\\n', 2333 maxlen=78, verb=1): 2334 """insert any '\\' chars for the given command 2335 - insert between start and end positions 2336 - apply specified wrap string wstring 2337 return a new string, in any case""" 2338 2339 global wrap_verb 2340 2341 if end < 0: end = len(command) - start - 1 2342 2343 nfirst = num_leading_line_spaces(command,start,1) # note initial indent 2344 prefix = get_next_indentation(command,start,end) 2345 sskip = nfirst # number of init spaces expected 2346 plen = len(prefix) 2347 newcmd = '' 2348 cur = start 2349 2350 if verb > 1: print("+d insert wrappers: nfirst=%d, prefix='%s', plen=%d" \ 2351 % (nfirst, prefix, plen)) 2352 2353 #pdb.set_trace() 2354 2355 # rewrite: create new command strings after each wrap 29 May 2009 2356 while needs_wrapper(command,maxlen,cur,end): 2357 endposn = command.find('\n',cur) 2358 if needs_wrapper(command,maxlen,cur,endposn): # no change on this line 2359 2360 lposn = find_last_space(command, cur+sskip, endposn, maxlen-sskip) 2361 2362 # if the last space is farther in than next indent, wrap 2363 # (adjust initial skip for any indent) 2364 if sskip+cur < lposn: # woohoo, wrap away (at lposn) 2365 newcmd = newcmd + command[cur:lposn+1] + wstring 2366 # modify command to add prefix, reset end and cur 2367 command = prefix + command[lposn+1:] 2368 end = end + plen - (lposn+1) 2369 sskip = nfirst + plen # now there is a prefix to skip 2370 cur = 0 2371 continue 2372 2373 # no change: 2374 # either the line does not need wrapping, or there is no space to do it 2375 if endposn < 0: endposn = end # there may not be a '\n' 2376 newcmd += command[cur:endposn+1] 2377 cur = endposn + 1 2378 2379 if cur <= end: newcmd += command[cur:end+1] # add remaining string 2380 2381 return newcmd 2382 2383def get_next_indentation(command,start=0,end=-1): 2384 """get any '#' plus leading spaces, from beginning or after first '\\\n'""" 2385 if end < 0: end = len(command) - start - 1 2386 2387 spaces = num_leading_line_spaces(command,start,1) 2388 prefix = command[start:start+spaces]+' ' # grab those spaces, plus 4 2389 # now check for an indention prefix 2390 posn = command.find('\\\n', start) 2391 pn = command.find('\n', start) # but don't continue past current line 2392 if posn >= 0 and posn < pn: 2393 spaces = num_leading_line_spaces(command,posn+2,1) 2394 if posn > start and spaces >= 2: 2395 prefix = command[posn+2:posn+2+spaces] # grab those spaces 2396 2397 return prefix 2398 2399def needs_wrapper(command, maxlen=78, start=0, end=-1): 2400 """does the current string need line wrappers 2401 2402 a string needs wrapping if there are more than 78 characters between 2403 any previous newline, and the next newline, wrap, or end""" 2404 2405 if end < 0: end_posn = len(command) - 1 2406 else: end_posn = end 2407 2408 cur_posn = start 2409 remain = end_posn - cur_posn 2410 while remain > maxlen: 2411 2412 # find next '\\\n' 2413 posn = command.find('\\\n', cur_posn) 2414 if 0 <= posn-cur_posn <= maxlen: # adjust and continue 2415 cur_posn = posn + 2 2416 remain = end_posn - cur_posn 2417 continue 2418 2419 # else find next '\n' 2420 posn = command.find('\n', cur_posn) 2421 if 0 <= posn-cur_posn <= maxlen: # adjust and continue 2422 cur_posn = posn + 1 2423 remain = end_posn - cur_posn 2424 continue 2425 2426 # otherwise, space means wrap, else not 2427 if find_next_space(command, cur_posn, 1) > cur_posn: return 1 2428 return 0 2429 2430 return 0 # if we get here, line wrapping is not needed 2431 2432def find_command_end(command, start=0): 2433 """find the next '\n' that is not preceeded by '\\', or return the 2434 last valid position (length-1)""" 2435 2436 length = len(command) 2437 end = start-1 2438 while 1: 2439 start = end + 1 2440 end = command.find('\n',start) 2441 2442 if end < 0: return length-1 # not in command 2443 elif end > start and command[end-1] == '\\': 2444 if length > end+1 and command[start] == '#' \ 2445 and command[end+1] != '#': 2446 return end # since comments cannot wrap 2447 else: continue 2448 return end # found 2449 2450def num_leading_line_spaces(istr,start,pound=0): 2451 """count the number of leading non-whitespace chars 2452 (newline chars are not be counted, as they end a line) 2453 if pound, skip any leading '#'""" 2454 2455 length = len(istr) 2456 if start < 0: start = 0 2457 if length < 1 or length <= start: return 0 2458 posn = start 2459 if pound and istr[posn] == '#': posn += 1 2460 2461 while posn < length and istr[posn].isspace() and istr[posn] != '\n': 2462 posn += 1 2463 2464 if posn == length: return 0 # none found 2465 return posn-start # index equals num spaces from start 2466 2467def find_next_space(istr,start,skip_prefix=0): 2468 """find (index of) first space after start that isn't a newline 2469 (skip any leading indendation if skip_prefix is set) 2470 return -1 if none are found""" 2471 2472 length = len(istr) 2473 index = start 2474 if skip_prefix: index += num_leading_line_spaces(istr,start,1) 2475 2476 while 1: 2477 if index >= length: break 2478 if istr[index] != '\n' and istr[index].isspace(): break 2479 index += 1 2480 2481 if index >= length : return -1 2482 return index 2483 2484def find_last_space(istr,start,end,max_len=-1,stretch=1): 2485 """find (index of) last space in current line range that isn't a newline 2486 if stretch and not found, search towards end 2487 return start-1 if none are found""" 2488 2489 if end < 0: end = len(istr) - 1 2490 if max_len >= 0 and end-start >= max_len: index = start+max_len-1 2491 else: index = end 2492 2493 posn = index # store current position in case of stretch 2494 2495 while posn >= start and (istr[posn] == '\n' or not istr[posn].isspace()): 2496 posn -= 1 2497 2498 if posn < start and stretch: # then search towards end 2499 posn = index 2500 while posn <= end and (istr[posn] == '\n' or not istr[posn].isspace()): 2501 posn += 1 2502 if posn > end: posn = start-1 # still failed 2503 2504 return posn # for either success or failure 2505 2506def nuke_final_whitespace(script, skipchars=[' ', '\t', '\n', '\\'], 2507 append='\n\n'): 2508 """replace final white characters with newlines""" 2509 slen = len(script) 2510 ind = slen-1 2511 while ind > 0 and script[ind] in skipchars: ind -= 1 2512 2513 return script[0:ind+1]+append 2514 2515# end line_wrapper functions 2516# ---------------------------------------------------------------------- 2517 2518# ---------------------------------------------------------------------- 2519# other functions 2520 2521# 17 May, 2008 [rickr] 2522def vals_are_multiples(num, vals, digits=4): 2523 """decide whether every value in 'vals' is a multiple of 'num' 2524 (vals can be a single float or a list of them) 2525 2526 Note, 'digits' can be used to specify the number of digits of accuracy 2527 in the test to see if a ratio is integral. For example: 2528 vals_are_multiples(1.1, 3.3001, 3) == 1 2529 vals_are_multiples(1.1, 3.3001, 4) == 0 2530 2531 return 1 if true, 0 otherwise (including error)""" 2532 2533 if num == 0.0: return 0 2534 2535 try: 2536 l = len(vals) 2537 vlist = vals 2538 except: 2539 vlist = [vals] 2540 2541 for val in vlist: 2542 rat = val/num 2543 rem = rat - int(rat) 2544 2545 if round(rem,digits) != 0.0: return 0 2546 2547 return 1 2548 2549def vals_are_constant(vlist, cval=None): 2550 """vals_are_constant(vlist, cval=None) 2551 2552 determine whether every value in vlist is equal to cval 2553 (if cval == None, use vlist[0])""" 2554 2555 if vlist == None: return 1 2556 if len(vlist) < 2: return 1 2557 2558 if cval == None: cval = vlist[0] 2559 2560 for val in vlist: 2561 if val != cval: return 0 2562 return 1 2563 2564def vals_are_positive(vlist): 2565 """determine whether every value in vlist is positive""" 2566 for val in vlist: 2567 if val <= 0: return 0 2568 return 1 2569 2570def vals_are_0_1(vlist): 2571 """determine whether every value in vlist is either 0 or 1""" 2572 for val in vlist: 2573 if val != 0 and val != 1: return 0 2574 return 1 2575 2576def vals_are_sorted(vlist, reverse=0): 2577 """determine whether values non-decreasing (or non-inc if reverse)""" 2578 if vlist == None: return 1 2579 if len(vlist) < 2: return 1 2580 2581 rval = 1 2582 try: 2583 for ind in range(len(vlist)-1): 2584 if reverse: 2585 if vlist[ind] < vlist[ind+1]: 2586 rval = 0 2587 break 2588 else: 2589 if vlist[ind] > vlist[ind+1]: 2590 rval = 0 2591 break 2592 except: 2593 print("** failed to detect sorting in list: %s" % vlist) 2594 rval = 0 2595 2596 return rval 2597 2598def vals_are_increasing(vlist, reverse=0): 2599 """determine whether values strictly increasing (or dec if reverse)""" 2600 if vlist == None: return 1 2601 if len(vlist) < 2: return 1 2602 2603 rval = 1 2604 try: 2605 for ind in range(len(vlist)-1): 2606 if reverse: 2607 if vlist[ind] <= vlist[ind+1]: 2608 rval = 0 2609 break 2610 else: # verify increasing 2611 if vlist[ind] >= vlist[ind+1]: 2612 rval = 0 2613 break 2614 except: 2615 print("** failed to detect sorting in list: %s" % vlist) 2616 rval = 0 2617 2618 return rval 2619 2620def vals_are_unique(vlist, dosort=1): 2621 """determine whether (possibly unsorted) values are unique 2622 - use memory to go for N*log(N) speed""" 2623 2624 if vlist == None: return 1 2625 if len(vlist) < 2: return 1 2626 2627 # copy and sort 2628 dupe = vlist[:] 2629 if dosort: dupe.sort() 2630 2631 rval = 1 2632 try: 2633 for ind in range(len(dupe)-1): 2634 if dupe[ind] == dupe[ind+1]: 2635 rval = 0 2636 break 2637 except: 2638 print("** uniq: failed to compare list elements in %s" % vlist) 2639 rval = 0 2640 2641 del(dupe) 2642 2643 return rval 2644 2645def lists_are_same(list1, list2, epsilon=0, doabs=0): 2646 """return 1 if the lists have similar values, else 0 2647 2648 similar means difference <= epsilon 2649 """ 2650 if not list1 and not list2: return 1 2651 if not list1: return 0 2652 if not list2: return 0 2653 if len(list1) != len(list2): return 0 2654 2655 for ind in range(len(list1)): 2656 if doabs: 2657 v1 = abs(list1[ind]) 2658 v2 = abs(list2[ind]) 2659 else: 2660 v1 = list1[ind] 2661 v2 = list2[ind] 2662 if v1 != v2: return 0 2663 if epsilon: 2664 if abs(v1-v2) > epsilon: return 0 2665 2666 return 1 2667 2668def list_intersect(listA, listB, sort=1): 2669 """return a list of the elements that are in both lists 2670 2671 if sort, sort the result 2672 """ 2673 2674 # if either is empty, the match list is as well 2675 if not listA or not listB: return [] 2676 2677 rlist = [v for v in listA if v in listB] 2678 2679 if sort: rlist.sort() 2680 2681 return rlist 2682 2683def list_diff(listA, listB, dtype='A-B', sort=1): 2684 """return a list of the elements differ between the lists 2685 2686 dtype: A-B : return elements in listA that are not list B 2687 B-A : return elements in listB that are not list A 2688 all : return all diffs (same as A-B and B-A) 2689 2690 return a list of newly created elements 2691 """ 2692 2693 # keep logic simple 2694 2695 if dtype == 'A-B': 2696 rlist = [v for v in listA if v not in listB] 2697 elif dtype == 'B-A': 2698 rlist = [v for v in listB if v not in listA] 2699 else: 2700 rlist = [v for v in listA if v not in listB] 2701 rlist.extend([v for v in listB if v not in listA]) 2702 2703 if sort: 2704 rlist.sort() 2705 2706 return rlist 2707 2708def string_to_float_list(fstring): 2709 """return a list of floats, converted from the string 2710 return None on error 2711 """ 2712 2713 if type(fstring) != str: return None 2714 slist = fstring.split() 2715 2716 if len(slist) == 0: return [] 2717 2718 try: flist = [float(sval) for sval in slist] 2719 except: return None 2720 2721 return flist 2722 2723def string_to_type_list(sdata, dtype=float): 2724 """return a list of dtype, converted from the string 2725 return None on error 2726 """ 2727 2728 if type(sdata) != str: return None 2729 slist = sdata.split() 2730 2731 if len(slist) == 0: return [] 2732 2733 # if going to int, use float as an intermediate step 2734 if dtype == int: 2735 try: slist = [float(sval) for sval in slist] 2736 except: return None 2737 2738 try: dlist = [dtype(sval) for sval in slist] 2739 except: return None 2740 2741 return dlist 2742 2743def float_list_string(vals, nchar=7, ndec=3, nspaces=2, mesg='', left=0): 2744 """return a string to display the floats: 2745 vals : the list of float values 2746 nchar : [7] number of characters to display per float 2747 ndec : [3] number of decimal places to print to 2748 nspaces : [2] number of spaces between each float 2749 """ 2750 2751 if left: form = '%-*.*f%*s' 2752 else: form = '%*.*f%*s' 2753 2754 istr = mesg 2755 for val in vals: istr += form % (nchar, ndec, val, nspaces, '') 2756 2757 return istr 2758 2759def gen_float_list_string(vals, mesg='', nchar=0, left=0): 2760 """mesg is printed first, if nchar>0, it is min char width""" 2761 2762 istr = mesg 2763 2764 if left: form = '%-' 2765 else: form = '%' 2766 2767 if nchar > 0: 2768 form += '*g ' 2769 for val in vals: istr += form % (nchar, val) 2770 else: 2771 form += 'g ' 2772 for val in vals: istr += form % val 2773 2774 return istr 2775 2776def int_list_string(ilist, mesg='', nchar=0, sepstr=' '): 2777 """like float list string, but use general printing 2778 (mesg is printed first, if nchar>0, it is min char width)""" 2779 2780 istr = mesg 2781 2782 if nchar > 0: slist = ['%*d' % (nchar, val) for val in ilist] 2783 else: slist = ['%d' % val for val in ilist] 2784 istr += sepstr.join(slist) 2785 2786 return istr 2787 2788def invert_int_list(ilist, top=-1, bot=0): 2789 """invert the integer list with respect to bot and top 2790 i.e. return a list of integers from bot to top that are not in 2791 the passed list 2792 """ 2793 if top < bot: 2794 print('** invert_int_list requires bot<=top (have %d, %d)' % (bot, top)) 2795 return [] 2796 2797 return [ind for ind in range(bot, top+1) if not ind in ilist] 2798 2799def is_valid_int_list(ldata, imin=0, imax=-1, whine=0): 2800 """check whether: 2801 o ldata is a of type [] 2802 o values are of type int 2803 o values are in within imin..imax (only if imin <= imax) 2804 o if whine: complain on error 2805 return 1 on true, 0 on false""" 2806 2807 if not ldata or type(ldata) != type([]): 2808 if whine: print("** not valid as a list: '%s'" % ldata) 2809 2810 for ind in range(len(ldata)): 2811 val = ldata[ind] 2812 if type(val) != type(0): 2813 if whine: print("** non-int value %d in int list (@ %d)" % (val,ind)) 2814 return 0 2815 if imin <= imax: # then also test bounds 2816 if val < imin: 2817 if whine: print("** list value %d not in [%d,%d]" %(val,imin,imax)) 2818 return 0 2819 elif val > imax: 2820 if whine: print("** list value %d not in [%d,%d]" %(val,imin,imax)) 2821 return 0 2822 return 1 2823 2824def data_to_hex_str(data): 2825 """convert raw data to hex string in groups of 4 bytes""" 2826 2827 if not data: return '' 2828 2829 dlen = len(data) # total length in bytes 2830 groups = (dlen+3) // 4 # number of 4-byte blocks to create 2831 remain = dlen 2832 retstr = '' # return string 2833 2834 for group in range(groups): 2835 if group > 0: retstr += ' ' 2836 retstr += '0x' 2837 if remain >= 4: llen = 4 2838 else: llen = remain 2839 2840 for ind in range(llen): 2841 retstr += '%02x' % data[dlen-remain+ind] 2842 2843 remain -= llen 2844 2845 return retstr 2846 2847def section_divider(hname='', maxlen=74, hchar='=', endchar=''): 2848 """return a title string of 'hchar's with the middle chars set to 'hname' 2849 if endchar is set, put at both ends of header 2850 e.g. section_divider('volreg', endchar='##') """ 2851 if len(hname) > 0: name = ' %s ' % hname 2852 else: name = '' 2853 2854 if endchar != '': maxlen -= 2*len(endchar) 2855 rmlen = len(name) 2856 if rmlen >= maxlen: 2857 print("** S_DIV_STR, rmlen=%d exceeds maxlen=%d" % (rmlen, maxlen)) 2858 return name 2859 prelen = (maxlen - rmlen) // 2 # basically half the chars 2860 postlen = maxlen - rmlen - prelen # other 'half' 2861 2862 return endchar + prelen*hchar + name + postlen*hchar + endchar 2863 2864def get_command_str(args=[], preamble=1, comment=1, quotize=1, wrap=1): 2865 """return a script generation command 2866 2867 args: arguments to apply 2868 preample: start with "script generated by..." 2869 comment: have text '#' commented out 2870 quotize: try to quotize any arguments that need it 2871 wrap: add line wrappers 2872 """ 2873 2874 if args == []: args = sys.argv 2875 2876 if preamble: hstr = '\n# %s\n# script generated by the command:\n#\n' \ 2877 % section_divider() 2878 else: hstr = '' 2879 2880 if comment: cpre = '# ' 2881 else: cpre = '' 2882 2883 # note command and args 2884 cmd = os.path.basename(args[0]) 2885 if quotize: args = ' '.join(quotize_list(args[1:],'')) 2886 else: args = ' '.join(args[1:],'') 2887 2888 cstr = '%s%s%s %s\n' % (hstr, cpre, cmd, args) 2889 2890 if wrap: return add_line_wrappers(cstr) 2891 else: return cstr 2892 2893def max_len_in_list(vlist): 2894 mval = 0 2895 try: 2896 for v in vlist: 2897 if len(v) > mval: mval = len(v) 2898 except: 2899 print('** max_len_in_list: cannot compute lengths') 2900 return mval 2901 2902def get_rank(data, style='dense', reverse=0, uniq=0): 2903 """return the rank order of indices given values, 2904 i.e. for each value, show its ordered index 2905 e.g. 3.4 -0.3 4.9 2.0 ==> 2 0 3 1 2906 2907 Sort the data, along with indices, then repeat on the newly 2908 sorted indices. If not uniq, set vals to minimum for repeated set. 2909 Note that sorting lists sorts on both elements, so first is smallest. 2910 2911 styles: (e.g. given input -3 5 5 6, result is ...) 2912 2913 competition: competition ranking (result 1 2 2 4 ) 2914 dense: dense ranking (result 1 2 2 3 ) 2915 {also default from 3dmerge and 3dRank} 2916 2917 return status (0=success) and the index order 2918 """ 2919 2920 dlen = len(data) 2921 2922 # maybe reverse the sort order 2923 if reverse: 2924 maxval = max(data) 2925 dd = [maxval-val for val in data] 2926 else: dd = data 2927 2928 # sort data and original position 2929 dd = [[dd[ind], ind] for ind in range(dlen)] 2930 dd.sort() 2931 2932 # invert postion list by repeating above, but with index list as data 2933 # (bring original data along for non-uniq case) 2934 dd = [[dd[ind][1], ind, dd[ind][0]] for ind in range(dlen)] 2935 2936 # deal with repeats: maybe modify d[1] from ind, depending on style 2937 if not uniq: 2938 if style == 'dense': 2939 cind = dd[0][1] # must be 0 2940 for ind in range(dlen-1): # compare next to current 2941 if dd[ind+1][2] == dd[ind][2]: 2942 dd[ind+1][1] = cind 2943 else: 2944 cind += 1 2945 dd[ind+1][1] = cind 2946 elif style == 'competition': 2947 for ind in range(dlen-1): # compare next to current 2948 if dd[ind+1][2] == dd[ind][2]: 2949 dd[ind+1][1] = dd[ind][1] 2950 else: 2951 print("** UTIL.GR: invalid style '%s'" % style) 2952 return 1, [] 2953 2954 dd.sort() 2955 2956 return 0, [dd[ind][1] for ind in range(dlen)] 2957 2958# ---------------------------------------------------------------------- 2959# wildcard construction functions 2960# ---------------------------------------------------------------------- 2961 2962def first_last_match_strs(slist): 2963 """given a list of strings, return the first and last consistent strings 2964 (i.e. all strings have the form first*last) 2965 2966 e.g. given ['subj_A1.txt', 'subj_B4.txt', 'subj_A2.txt' ] 2967 return 'subj_' and '.txt' 2968 """ 2969 2970 if type(slist) != list: 2971 print('** FL match strings requires a list') 2972 return '', '' 2973 2974 if not slist: return '', '' 2975 2976 maxlen = len(slist[0]) 2977 hmatch = maxlen # let them shrink 2978 tmatch = maxlen 2979 for sind in range(1, len(slist)): 2980 if slist[0] == slist[sind]: continue 2981 2982 hmatch = min(hmatch, len(slist[sind])) 2983 tmatch = min(tmatch, len(slist[sind])) 2984 2985 # find first left diff 2986 i = 0 2987 while i < hmatch: 2988 if slist[sind][i] != slist[0][i]: break 2989 i += 1 2990 hmatch = min(hmatch, i) 2991 2992 # find first right diff (index from 1) 2993 i = 1 2994 while i <= tmatch: 2995 if slist[sind][-i] != slist[0][-i]: break 2996 i += 1 2997 tmatch = min(tmatch, i-1) 2998 2999 if hmatch+tmatch > maxlen: # weird, but constructable 3000 tmatch = maxlen - hmatch # so shrink to fit 3001 3002 # list[-0:] is not empty but is the whole list 3003 if tmatch > 0: tstr = slist[0][-tmatch:] 3004 else: tstr = '' 3005 3006 return slist[0][0:hmatch], tstr 3007 3008def list_files_by_glob(L, sort=False, exit_on_miss=False) : 3009 """Input a list L of one or more strings to glob (fiiine, the input L 3010 can also be a string, which will just be made into a list 3011 immediately, L = [L]), and then glob each one, sort that piece's 3012 results, and append it to the list of results. Return the full 3013 list of files. 3014 3015 The program can be made to exist if any piece fails to find a 3016 result, by setting 'exit_on_miss=True'. 3017 3018 Sorting of the full/final list is NOT done by default (can be 3019 turned on with 'sort=True'). 3020 3021 """ 3022 3023 if not(L): return [] 3024 3025 if type(L) == str : L = [L] 3026 elif type(L) != list : BASE.EP("Input a list (or str)") 3027 3028 N = len(L) 3029 cbad = 0 3030 out = [] 3031 3032 for ii in range(N): 3033 sublist = L[ii].split() 3034 for item in sublist: 3035 gout = glob.glob(item) 3036 if not(gout) : 3037 cbad+= 1 3038 BASE.WP("No files for this part of list: {}".format(item)) 3039 else: 3040 gout.sort() 3041 out.extend(gout) 3042 3043 if cbad and exit_on_miss : 3044 BASE.EP("No findings for {} parts of this list. Bye.".format(cbad)) 3045 3046 if sort : 3047 out.sort() 3048 3049 return out 3050 3051def glob2stdout(globlist): 3052 """given a list of glob forms, print all matches to stdout 3053 3054 This is meant to be a stream workaround to shell errors 3055 like, "Argument list too long". 3056 3057 echo 'd1/*.dcm' 'd2/*.dcm' \\ 3058 | afni_python_wrapper.py -listfunc glob2stdout - 3059 """ 3060 for gform in globlist: 3061 for fname in glob.glob(gform): 3062 print(fname) 3063 3064def glob_form_from_list(slist): 3065 """given a list of strings, return a glob form 3066 3067 e.g. given ['subjA1.txt', 'subjB4.txt', 'subjA2.txt' ] 3068 return 'subj*.txt' 3069 3070 Somewhat opposite list_minus_glob_form(). 3071 """ 3072 3073 if len(slist) == 0: return '' 3074 if vals_are_constant(slist): return slist[0] 3075 3076 first, last = first_last_match_strs(slist) 3077 if not first and not last: return '' # failure 3078 globstr = '%s*%s' % (first,last) 3079 3080 return globstr 3081 3082def glob_form_matches_list(slist, ordered=1): 3083 """given a list of strings, make a glob form, and then test that against 3084 the actual files on disk 3085 3086 if ordered: files must match exactly (i.e. slist must be sorted) 3087 else: slist does not need to be sorted 3088 """ 3089 3090 slen = len(slist) 3091 3092 # check trivial cases of lengths 0 and 1 3093 if slen == 0: return 1 3094 if slen == 1: 3095 if os.path.isfile(slist[0]): return 1 3096 else: return 0 3097 3098 globstr = glob_form_from_list(slist) 3099 glist = glob.glob(globstr) 3100 glist.sort() 3101 3102 # quick check: lengths must match 3103 if len(glist) != slen: return 0 3104 3105 if ordered: 3106 inlist = slist 3107 else: 3108 inlist = slist[:] 3109 inlist.sort() 3110 3111 # now files must match exactly (between inlist and glist) 3112 for ind in range(slen): 3113 if glist[ind] != inlist[ind]: return 0 3114 3115 # they must match 3116 return 1 3117 3118 3119def list_minus_glob_form(inlist, hpad=0, tpad=0, keep_dent_pre=0, strip=''): 3120 """given a list of strings, return the inner part of the list that varies 3121 (i.e. remove the consistent head and tail elements) 3122 3123 e.g. given ['subjA1.txt', 'subjB4.txt', 'subjA2.txt' ] 3124 return [ 'A1', 'B4', 'A2' ] 3125 3126 hpad NPAD : number of characters to pad at prefix 3127 tpad NPAD : number of characters to pad at suffix 3128 keep_dent_pre Y/N : (flag) keep entire prefix from directory entry 3129 strip : one of ['', 'dir', 'file', 'ext', 'fext'] 3130 3131 If hpad > 0, then pad with that many characters back into the head 3132 element. Similarly, tpad pads forward into the tail. 3133 3134 e.g. given ['subjA1.txt', 'subjB4.txt', 'subjA2.txt' ] 3135 if hpad = 926 (or 4 :) and tpad = 1, 3136 return [ 'subjA1.', 'subjB4.', 'subjA2.' ] 3137 3138 If keep_dent_pre is set, then (if '/' is found) decrement hlen until 3139 that '/'. 3140 3141 e.g. given ['dir/subjA1.txt', 'dir/subjB4.txt', 'dir/subjA2.txt' ] 3142 -> return = [ 'A1.', 'B4.', 'A2.' ] 3143 with keep_dent_pre == 1: 3144 -> return = [ 'subjA1.', 'subjB4.', 'subjA2.' ] 3145 3146 Somewhat opposite glob_form_from_list(). 3147 """ 3148 3149 if len(inlist) <= 1: return inlist 3150 3151 # init with original 3152 slist = inlist 3153 3154 # maybe make a new list of stripped elements 3155 stripnames = ['dir', 'file', 'ext', 'fext'] 3156 if strip != '' and strip not in stripnames: 3157 print('** LMGF: bad strip %s' % strip) 3158 strip = '' 3159 3160 if strip in stripnames: 3161 ss = [] 3162 for inname in inlist: 3163 if strip == 'dir': 3164 dname, fname = os.path.split(inname) 3165 ss.append(fname) 3166 elif strip == 'file': 3167 dname, fname = os.path.split(inname) 3168 ss.append(dname) 3169 elif strip == 'ext': 3170 fff, ext = os.path.splittext(inname) 3171 ss.append(fff) 3172 elif strip == 'fext': 3173 fff, ext = os.path.splittext(inname) 3174 ss.append(fff) 3175 else: 3176 print('** LMGF: doubly bad strip %s' % strip) 3177 break 3178 # check for success 3179 if len(ss) == len(slist): slist = ss 3180 3181 if hpad < 0 or tpad < 0: 3182 print('** list_minus_glob_form: hpad/tpad must be non-negative') 3183 hpad = 0 ; tpad = 0 3184 3185 # get head, tail and note lengths 3186 head, tail = first_last_match_strs(slist) 3187 hlen = len(head) 3188 tlen = len(tail) 3189 3190 # adjust by padding, but do not go negative 3191 if hpad >= hlen: hlen = 0 3192 else: hlen -= hpad 3193 if tpad >= tlen: tlen = 0 3194 else: tlen -= tpad 3195 3196 # apply directory entry prefix, if requested 3197 if keep_dent_pre: 3198 s = slist[0] 3199 posn = s.rfind('/', 0, hlen) 3200 # if found, start at position to right of it 3201 # otherwise, use entire prefix 3202 if posn >= 0: hlen = posn + 1 3203 else: hlen = 0 3204 3205 # and return the list of center strings 3206 if tlen == 0: return [ s[hlen:] for s in slist ] 3207 else: return [ s[hlen:-tlen] for s in slist ] 3208 3209def glob_list_minus_pref_suf(pref, suf): 3210 """just strip the prefix and suffix from string list elements 3211 (for now, assume they are there) 3212 """ 3213 glist = glob.glob('%s*%s' % (pref, suf)) 3214 3215 plen = len(pref) 3216 slen = len(suf) 3217 3218 return [d[plen:-slen] for d in glist] 3219 3220def list_minus_pref_suf(slist, pref, suf, stripdir=1): 3221 """just strip the prefix and suffix from string list elements 3222 3223 if stripdir, remove leading directories 3224 3225 return status, stripped list 3226 3227 status = 0 : all strings have prefix and suffix 3228 1 : not all do 3229 -1 : on error 3230 """ 3231 3232 plen = len(pref) 3233 slen = len(suf) 3234 3235 # possibly strip of directory names 3236 if stripdir: 3237 flist = [] 3238 for sname in slist: 3239 dd, ff = os.path.split(sname) 3240 flist.append(ff) 3241 else: flist = slist 3242 3243 3244 rv = 0 3245 rlist = [] 3246 for fname in flist: 3247 if fname.startswith(pref): poff = plen 3248 else: poff = 0 3249 3250 if fname.endswith(suf): soff = slen 3251 else: soff = 0 3252 3253 if soff: rlist.append(fname[poff:-soff]) 3254 else: rlist.append(fname[poff:]) 3255 3256 if not poff or not soff: rv = 1 3257 3258 return rv, rlist 3259 3260def okay_as_lr_spec_names(fnames, verb=0): 3261 """check that names are okay as surface spec files, e.g. for afni_proc.py 3262 - must be 1 or 2 file names 3263 - must contain lh and rh, respectively 3264 - must otherwise be identical 3265 3266 if verb, output why failure occurs 3267 return 1 if okay, 0 otherwise 3268 """ 3269 nfiles = len(fnames) 3270 if nfiles == 0: return 1 # no problems, anyway 3271 if nfiles > 2: 3272 if verb: print('** only 1 or 2 spec files allowed (have %d)' % nfiles) 3273 return 0 3274 3275 if nfiles == 1: 3276 if fnames[0].find('lh')>=0 or fnames[0].find('rh')>=0: return 1 # success 3277 # failure 3278 if verb: print("** spec file '%s' missing 'lh' or 'rh'" % fnames[0]) 3279 return 0 3280 3281 # so we have 2 files 3282 3283 hlist = list_minus_glob_form(fnames, tpad=1) # go after following 'h' 3284 3285 for h in hlist: 3286 if h != 'rh' and h != 'lh': 3287 if verb: print('** multiple spec files must only differ in lh vs. rh') 3288 return 0 3289 3290 return 1 3291 3292def make_spec_var(fnames, vname='hemi'): 3293 """return a spec file variable and a list of replaced hemispheres 3294 3295 e.g. make_spec_var(['surface.lh.spec', 'surface.rh.spec']) returns 3296 surface.${hemi}.spec, ['lh', 'rh'] 3297 given 1 or 2 spec file names, return a single variable that 3298 represents them with $vname replacing the lh or rh 3299 3300 return '' on failure 3301 """ 3302 if not okay_as_lr_spec_names(fnames): return '', [] 3303 3304 nfiles = len(fnames) 3305 if nfiles == 0 or nfiles > 2: return '', [] 3306 3307 sfile = fnames[0] 3308 3309 if nfiles == 1: 3310 # just find lh or rh and replace it 3311 hh = 'lh' 3312 posn = sfile.find(hh) 3313 if posn < 0: 3314 hh = 'rh' 3315 posn = sfile.find(hh) 3316 if posn < 0: return '', [] # should not happen 3317 3318 return sfile[0:posn] + '${%s}'%vname + sfile[posn+2:], [hh] 3319 3320 # so nfiles == 2, use glob 3321 3322 head, tail = first_last_match_strs(fnames) 3323 hlen = len(head) 3324 3325 hemi = sfile[hlen:hlen+2] 3326 if hemi != 'lh' and hemi != 'rh': 3327 print('** MSV: bad lh/rh search from spec files: %s' % fnames) 3328 return '', [] 3329 3330 return sfile[0:hlen] + '${%s}'%vname + sfile[hlen+2:], ['lh', 'rh'] 3331 3332def parse_as_stim_list(flist): 3333 """parse filename list as PREFIX.INDEX.LABEL.SUFFIX, where the separators 3334 can be '.', '_' or nothing (though ignore PREFIX and SUFFIX, as well 3335 as those separators) 3336 3337 - strip PREFIX and SUFFIX (those are garbage) 3338 (if SUFFIX has a '.' after position 0, adjust the SUFFIX) 3339 - strip any leading digits as INDEXes 3340 3341 return Nx2 table where if one column entry is filled, they all are 3342 (None on failure for form a complete table) 3343 (note: blank INDEX or LABEL is okay, but they must all be) 3344 """ 3345 3346 if len(flist) < 1: return [] 3347 3348 # first get PREFIX and SUFFIX 3349 prefix, suffix = first_last_match_strs(flist) 3350 3351 # if suffix contains an extension, make the suffix into the extension 3352 dot = suffix.find('.') 3353 if dot < 0: dot = 0 3354 3355 # strip prefix, suffix: might include part of 'suffix' in label 3356 inner_list = list_minus_glob_form(flist, tpad=dot) 3357 3358 # then make table of the form <NUMBER><SEP><LABEL> 3359 s_table = [list(_parse_leading_int(name)) for name in inner_list] 3360 3361 # if any number does not exist, just return inner_list as LABELs 3362 for entry in s_table: 3363 if entry[0] < 0: return [[-1, label] for label in inner_list] 3364 3365 # return INDEX and LABEL (no SEP) 3366 return [[entry[0], entry[2]] for entry in s_table] 3367 3368def _parse_leading_int(name, seplist=['.','_','-']): 3369 """assuming name is a string starting with digits, parse name into 3370 val, sep, suffix 3371 3372 val = -1 if name does not start with a digit 3373 sep = one of {'.', '_', '-', ''} 3374 suffix = whatever remains after 'sep' 3375 """ 3376 3377 nlen = len(name) 3378 3379 if nlen < 1: return -1, '', '' 3380 3381 # first strip of any leading (non-negative) integer 3382 posn = 0 # count leading digits 3383 for char in name: 3384 if char.isdigit(): posn += 1 3385 else: break 3386 3387 if posn == 0: val = -1 3388 else: 3389 try: val = int(name[0:posn]) 3390 except: 3391 print("** _parse_leading_int: can't parse int from %s" % name) 3392 return 3393 3394 # if only a number, we're outta here 3395 if posn == nlen: return val, '', '' 3396 3397 # note any separator 3398 if name[posn] in seplist: 3399 sep = name[posn] 3400 posn += 1 3401 else: 3402 sep = '' 3403 3404 # aaaaaand, we're done 3405 return val, sep, name[posn:] 3406 3407def glob_form_has_match(form): 3408 """see if anything at all exists according to this glob form""" 3409 glist = glob.glob(form) 3410 glen = len(glist) 3411 del(glist) 3412 if glen > 0: return 1 3413 return 0 3414 3415def executable_dir(ename=''): 3416 """return the directory whre the ename program is located 3417 (by default, use argv[0])""" 3418 if ename == '': ee = sys.argv[0] 3419 else: ee = ename 3420 3421 dname = os.path.dirname(ee) 3422 return os.path.abspath(dname) 3423 3424def common_dir(flist): 3425 """return the directory name that is common to all files (unless trivial)""" 3426 dname, junk = first_last_match_strs(flist) 3427 if len(dname) > 0 and dname[-1] == '/': dname = dname[0:-1] 3428 if not os.path.isdir(dname): dname = os.path.dirname(dname) 3429 3430 if is_trivial_dir(dname): return '' 3431 return dname 3432 3433def common_parent_dirs(flists): 3434 """return parent directories 3435 3436 flists = lists of file names (each element should be a list) 3437 3438 return: 3439 top_dir (common to all parents (files), '' if not used) 3440 parent_dir (for each flist, common parent) 3441 short_dir (for each flist, common parent under top_dir) 3442 short_name (for each flist, file names under parent dirs) 3443 3444 if top_dir has at least 2 levels, use it 3445 """ 3446 if type(flists) != list: 3447 print('** common_parent_dirs: bad flists type') 3448 return None, None, None, None 3449 for ind in range(len(flists)): 3450 flist = flists[ind] 3451 if type(flist) != list: 3452 print('** common_parent_dirs: bad flist[%d] type' % ind) 3453 return None, None, None, None 3454 3455 # get top_dir and parents 3456 all_pars = [] 3457 par_dirs = [] 3458 short_names = [] 3459 for flist in flists: 3460 # track parent dirs 3461 parent = common_dir(flist) 3462 if parent == '/' or is_trivial_dir(parent): parent = '' 3463 par_dirs.append(parent) 3464 3465 # and make short names 3466 plen = len(parent) 3467 if plen > 0: start = plen+1 3468 else: start = 0 3469 short_names.append([fname[start:] for fname in flist]) 3470 3471 # top is common to all parents 3472 top_dir = common_dir(par_dirs) 3473 if top_dir.count('/') <= 1: top_dir = '' 3474 3475 # now get all short dir names, under top dir 3476 if top_dir == '': short_dirs = par_dirs 3477 else: short_dirs = [child_dir_name(top_dir, pdir) for pdir in par_dirs] 3478 3479 return top_dir, par_dirs, short_dirs, short_names 3480 3481def child_dir_name(parent, child): 3482 """return the child directory name truncated under the parent""" 3483 if parent == '' or child == '': return child 3484 plen = len(parent) 3485 clen = len(child) 3486 3487 if child[0:plen] != parent: return child # not a proper child 3488 3489 # return everything after separator 3490 if clen < plen + 2: return '.' # trivial as child 3491 else: return child[plen+1:] # remove parent portion 3492 3493def is_trivial_dir(dname): 3494 """input a string 3495 return 1 if dname is empty or '.' 3496 else return 0 3497 """ 3498 if dname == None: return 1 3499 if dname == '' or dname == '.' or dname == './' : return 1 3500 3501 return 0 3502 3503def flist_to_table_pieces(flist): 3504 """dissect a file list 3505 input: list of file names 3506 output: 3507 - common directory name 3508 - short name list (names after common directory) 3509 - glob string from short names 3510 note: short names will be new data, never just references to input 3511 """ 3512 if len(flist) == 0: return '', [], '' 3513 3514 ddir = common_dir(flist) 3515 dirlen = len(ddir) 3516 if dirlen > 0: snames = [dset[dirlen+1:] for dset in flist] 3517 else: snames = [dset[:] for dset in flist] 3518 3519 globstr = glob_form_from_list(snames) 3520 3521 return ddir, snames, globstr 3522 3523def get_ids_from_dsets(dsets, prefix='', suffix='', hpad=0, tpad=0, verb=1): 3524 """return a list of subject IDs corresponding to the datasets 3525 3526 Make sure we have afni_name objects to take the prefix from. 3527 If simple strings, turn into afni_names. 3528 3529 Try list_minus_glob_form on the datasets. If that fails, try 3530 on the directories. 3531 3532 prefix, suffix: attach these to the resulting IDs 3533 hpad, tpad: padding to pass to list_minus_glob_form 3534 3535 return None on failure 3536 """ 3537 if hpad < 0 or tpad < 0: 3538 print('** get_ids_from_dsets: will not apply negative padding') 3539 hpad, tpad = 0, 0 3540 3541 if len(dsets) == 0: return None 3542 3543 # be more aggressive, use dataset prefix names 3544 # dlist = [dset.split('/')[-1] for dset in dsets] 3545 if type(dsets[0]) == str: 3546 nlist = [BASE.afni_name(dset) for dset in dsets] 3547 elif isinstance(dsets[0], BASE.afni_name): 3548 nlist = dsets 3549 else: 3550 print('** GIFD: invalid type for dset list, have value %s' % dsets[0]) 3551 return None 3552 3553 dlist = [dname.prefix for dname in nlist] 3554 3555 # if nothing to come from file prefixes, try the complete path names 3556 if vals_are_constant(dlist): dlist = dsets 3557 3558 slist = list_minus_glob_form(dlist, hpad, tpad) 3559 3560 # do some error checking 3561 for val in slist: 3562 if '/' in val: # no directories 3563 if verb > 0: print('** GIFD: IDs would have directories') 3564 return None 3565 3566 if len(slist) != len(dsets): # appropriate number of entries 3567 if verb > 0: print('** GIFD: length mis-match getting IDs') 3568 return None 3569 3570 if not vals_are_unique(slist): 3571 if verb > 0: print('** GIFD: final IDs are not unique') 3572 return None 3573 3574 return slist 3575 3576def insensitive_word_pattern(word): 3577 """replace each alphabetical char with [Ul], an upper/lower search pair 3578 use of 'either' suggested by G Irving at stackoverflow.com 3579 """ 3580 def either(c): 3581 if c.isalpha: return '[%s%s]'%(c.lower(),c.upper()) 3582 else: return c 3583 return ''.join(map(either,word)) 3584 3585def insensitive_glob(pattern): 3586 """return glob.glob, but where every alphabetic character is 3587 replaced by lower/upper pair tests 3588 """ 3589 return glob.glob(insensitive_word_pattern(pattern)) 3590 3591 3592def search_path_dirs(word, mtype=0, casematch=1): 3593 """return status and list of matching files 3594 3595 Could be more efficient, esp. with mtype=exact and casematch set, but 3596 I will strive for simplicity and consistency and see how it goes. 3597 3598 mtype : 0 = match any sub-word (i.e. look for DIR/*word*) 3599 1 = exact match (i.e. no wildcard, look for DIR/word) 3600 2 = prefix match (i.e. look for DIR/word*) 3601 casematch : flag: if set, case must match 3602 else, 'word' letters can be either case 3603 """ 3604 try: 3605 plist = os.environ['PATH'].split(os.pathsep) 3606 except: 3607 print('** search_path_dirs: no PATH var') 3608 return 1, [] 3609 3610 # if no casematch, look for upper/lower pairs 3611 if casematch: wpat = word 3612 else: wpat = insensitive_word_pattern(word) 3613 3614 # if not exact, surround with wildcard pattern 3615 if mtype == 0: form = '%s/*%s*' # any sub-word 3616 elif mtype == 1: form = '%s/%s' # exact match 3617 elif mtype == 2: form = '%s/%s*' # prefix match 3618 else: 3619 print('** search_path_dirs: illegal mtype = %s' % mtype) 3620 return 1, [] 3621 3622 # now just search for matches 3623 rlist = [] 3624 for pdir in plist: 3625 glist = glob.glob(form % (pdir, wpat)) 3626 glist.sort() 3627 if len(glist) > 0: rlist.extend(glist) 3628 3629 # make a new list based on os.path.realpath, to avoid links 3630 rlist = [os.path.realpath(pfile) for pfile in rlist] 3631 3632 return 0, get_unique_sublist(rlist) 3633 3634def which(pname): 3635 """like unix which command: return the first 'pname' in PATH 3636 (this is a simplified version of search_path_dirs()) 3637 3638 return a simple string, empty on error 3639 """ 3640 try: 3641 plist = os.environ['PATH'].split(os.pathsep) 3642 except: 3643 print('** which_prog: no PATH var') 3644 return '' 3645 3646 for pdir in plist: 3647 # accept pname having a path 3648 search = os.path.join(pdir, pname) 3649 if os.path.isfile(search) and os.access(search, os.X_OK): 3650 return search 3651 3652 return '' 3653 3654 3655def num_found_in_path(word, mtype=0, casematch=1): 3656 """a simple wrapper to print search_path_dirs results 3657 3658 Search for given 'word' in path, and print out list elements 3659 with element prefix of 'indent'. 3660 3661 mtype : 0 = match any sub-word (i.e. look for DIR/*word*) 3662 1 = exact match (i.e. no wildcard, look for DIR/word) 3663 2 = prefix match (i.e. look for DIR/word*) 3664 casematch : flag: if set, case must match 3665 else, 'word' letters can be either case 3666 indent : prefix/separator string for list elements 3667 """ 3668 rv, rlist = search_path_dirs(word, mtype=mtype, casematch=casematch) 3669 if rv: return 0 3670 return len(rlist) 3671 3672def show_found_in_path(word, mtype=0, casematch=1, indent='\n '): 3673 """a simple wrapper to print search_path_dirs results 3674 3675 Search for given 'word' in path, and print out list elements 3676 with element prefix of 'indent'. 3677 3678 mtype : 0 = match any sub-word (i.e. look for DIR/*word*) 3679 1 = exact match (i.e. no wildcard, look for DIR/word) 3680 2 = prefix match (i.e. look for DIR/word*) 3681 casematch : flag: if set, case must match 3682 else, 'word' letters can be either case 3683 indent : prefix/separator string for list elements 3684 """ 3685 rv, rlist = search_path_dirs(word, mtype=mtype, casematch=casematch) 3686 if not rv: print(indent+indent.join(rlist)) 3687 3688# ---------------------------------------------------------------------- 3689# mathematical functions: 3690# vector routines: sum, sum squares, mean, demean 3691# ---------------------------------------------------------------------- 3692 3693def loc_sum(vals): 3694 """in case 'sum' does not exist, such as on old machines""" 3695 3696 try: tot = sum(vals) 3697 except: 3698 tot = 0 3699 for val in vals: tot += val 3700 return tot 3701 3702def sumsq(vals): 3703 """return the sum of the squared values""" 3704 ssq = 0 3705 for val in vals: ssq += (val*val) 3706 return ssq 3707 3708def euclidean_norm(vals): 3709 """name is toooooo long""" 3710 3711 if len(vals) < 1: return 0.0 3712 return math.sqrt(sumsq(vals)) 3713 3714def L2_norm(vals): 3715 return euclidean_norm(vals) 3716 3717def weighted_enorm(vals, weights): 3718 3719 if len(vals) < 1: return 0.0 3720 if len(vals) != len(weights): return 0.0 3721 3722 sum = 0.0 3723 for ind in range(len(vals)): 3724 vv = vals[ind]*weights[ind] 3725 sum += vv*vv 3726 return math.sqrt(sum) 3727 3728def dotprod(v1,v2): 3729 """compute the dot product of 2 vectors""" 3730 try: dsum = loc_sum([v1[i]*v2[i] for i in range(len(v1))]) 3731 except: 3732 print('** cannot take dotprod() of these elements') 3733 dsum = 0 3734 return dsum 3735 3736def affine_to_params_6(avec, verb=1): 3737 """convert rotation/shift affine "matrix" to 6 parameters 3738 (e.g. as in 3dvolreg 1Dmatrix format to 1Dfile format) 3739 3740 matvec: length 12+ vector (row major order) 3741 return: length 6 param vector: 3742 roll, pitch, yaw, dx, dy, dz 3743 """ 3744 3745 rvec = [0.0]*6 3746 3747 if len(avec) < 12: 3748 print('** affine_to_params_6: requires length 12+ vector, have %d' \ 3749 % len(avec)) 3750 return rvec 3751 3752 # rotations 3753 rvec[0] = 180.0/math.pi * math.atan2(avec[9], avec[10]) 3754 rvec[1] = 180.0/math.pi *-math.asin (avec[8]) 3755 rvec[2] = 180.0/math.pi * math.atan2(avec[4], avec[0]) 3756 3757 # deltas 3758 rvec[3] = avec[3] 3759 rvec[4] = avec[7] 3760 rvec[5] = avec[11] 3761 3762 return rvec 3763 3764def maxabs(vals): 3765 """convenience function for the maximum of the absolute values""" 3766 if len(vals) == 0: return 0 3767 return max([abs(v) for v in vals]) 3768 3769def ndigits_lod(num, base=10): 3770 """return the number of digits to the left of the decimal""" 3771 anum = abs(num) 3772 if base == 10: return 1+int(math.log10(anum)) 3773 else: return 1+int(math.log10(anum)/math.log10(base)) 3774 3775# almost identical to demean, but just return the mean 3776def mean(vec, ibot=-1, itop=-1): 3777 """return the vector mean, from index ibot to index itop 3778 3779 if ibot == -1, ibot will be 0 3780 if itop == -1, itop will be len-1""" 3781 3782 if not vec: return 0.0 3783 if ibot > itop: 3784 print('** afni_util.mean: ibot (%d) > itop (%d)' % (ibot, itop)) 3785 return 0.0 3786 3787 vlen = len(vec) 3788 3789 if ibot < 0: ibot = 0 3790 if ibot > vlen-1: ibot = vlen-1 3791 if itop < 0: itop = vlen-1 3792 if itop > vlen-1: itop = vlen-1 3793 3794 tot = 0.0 3795 for ind in range(ibot,itop+1): 3796 tot += vec[ind] 3797 3798 return tot/(itop-ibot+1) 3799 3800 3801# convert from degrees to chord length 3802def deg2chordlen(theta, radius=1.0): 3803 """deg2chord_length(theta, radius=1.0): 3804 3805 Given theta in degrees (0<=theta<=180) and a radius (>=0), return the 3806 length of the chord with an arc that subtends central angle theta. 3807 (So the chord corresponds with an arc of a circle, and the arc subtends 3808 central angle theta.) 3809 This might help estimate motion distance due to a rotation. 3810 3811 For a circle of radius R and a central angle T in degrees, compute the 3812 resulting chord length (distance between points on the circle for the 3813 corresponding arc). If the chord has endpoints A and B, we are looking 3814 for the length of the segment (AB). 3815 3816 Note that if a perpendicular is dropped from the circle's center to AB, 3817 cutting it in half (call the length h), then we have: 3818 3819 sin(T/2) = h/R, so h = R * sin(T/2) 3820 3821 return 2*h (to get the entire chord length) 3822 """ 3823 3824 # put theta in [0,180] 3825 if theta < 0.0: theta = abs(theta) 3826 if theta > 360: theta = theta % 360 3827 if theta > 180: theta = 180 - theta 3828 3829 # ignore a negative radius 3830 if radius <= 0.0: return 0.0 3831 3832 # math.tan takes input in radians 3833 theta = theta * math.pi / 180.0 3834 3835 return 2.0 * radius * math.sin(theta/2.0) 3836 3837# ---------------------------------------------------------------------- 3838# vector manipulation functions 3839# ---------------------------------------------------------------------- 3840 3841# almost identical to mean, but subtract the mean instead of returning it 3842def demean(vec, ibot=-1, itop=-1): 3843 """demean the vector (in place), from index ibot to index itop 3844 3845 if ibot == -1, ibot will be 0 3846 if itop == -1, itop will be len-1 3847 3848 return 0 on success, 1 on error""" 3849 3850 if not vec: return 0 3851 if ibot > itop: 3852 print('** afni_util.demean: ibot (%d) > itop (%d)' % (ibot, itop)) 3853 return 1 3854 3855 vlen = len(vec) 3856 3857 if ibot < 0: ibot = 0 3858 if ibot > vlen-1: ibot = vlen-1 3859 if itop < 0: itop = vlen-1 3860 if itop > vlen-1: itop = vlen-1 3861 3862 # first compute the mean 3863 tot = 0.0 3864 for ind in range(ibot,itop+1): 3865 tot += vec[ind] 3866 mm = tot/(itop-ibot+1) 3867 3868 # now subract it 3869 for ind in range(ibot,itop+1): 3870 vec[ind] -= mm 3871 3872 return vec 3873 3874def lin_vec_sum(s1, vec1, s2, vec2): 3875 """return s1*[vec1] + s2*[vec2] 3876 note: vec2 can be None""" 3877 3878 if vec2 == None: 3879 return [s1*vec1[i] for i in range(len(vec1))] 3880 3881 l1 = len(vec1) 3882 l2 = len(vec2) 3883 if l1 != l2: 3884 print('** LVC: vectors have different lengths (%d, %d)' % (l1, l2)) 3885 return [] 3886 3887 return [s1*vec1[i]+s2*vec2[i] for i in range(l1)] 3888 3889def proj_onto_vec(v1, v2, unit_v2=0): 3890 """return vector v1 projected onto v2 3891 3892 unit_v2: flag denoting whether v2 is a unit vector 3893 3894 return <v1,v2>/<v2,v2> * v2""" 3895 3896 if unit_v2: scalar = dotprod(v1,v2) 3897 else: 3898 len2 = L2_norm(v2,v2) 3899 if len2 == 0: 3900 print('** cannot project onto <0> vector') 3901 return [] 3902 scalar = dotprod(v1,v2) * 1.0 / len2 3903 3904 return lin_vec_sum(scalar, v2, 0, None) 3905 3906def proj_out_vec(v1, v2, unit_v2=0): 3907 """return vector v1 with v2 projected out 3908 (return the component of v1 that is orthogonal to v2) 3909 3910 unit_v2: flag denoting whether v2 is a unit vector 3911 3912 return v1 - (v1 projected onto v2) 3913 3914 Note: y - proj(y,p) 3915 = y - <y,p>/<p,p> * pT = y - yTp/pTp * pT 3916 = y - <y,p>/|p|^2 * pT 3917 = y - <y,p>*(pTp)^-1 * pT (where (pTp)^-1*pT = pseudo-inverse) 3918 = (I - p (pTp)^-1 * pT) * y 3919 """ 3920 3921 return lin_vec_sum(1, v1, -1, proj_onto_vec(v1, v2, unit_v2)) 3922 3923# ---------------------------------------------------------------------- 3924# statistical routines - stdev, variance, ttest 3925# ---------------------------------------------------------------------- 3926 3927def stat_mean_abs_dev(data): 3928 """return the mean absolute deviation""" 3929 3930 if not data: return 0 3931 length = len(data) 3932 # length == 1 has MAD 0 3933 if length <= 1: return 0 3934 3935 if type(data[0]) == str: 3936 try: dd = [float(val) for val in data] 3937 except: 3938 print('** bad data for min_mean_max_stdev') 3939 return 0, 0, 0, 0 3940 else: dd = data 3941 3942 meanval = loc_sum(dd)/float(length) 3943 dsum = 0.0 3944 for val in dd: 3945 dsum += abs(val-meanval) 3946 return 1.0*dsum/length 3947 3948def min_mean_max_stdev(data): 3949 """return 4 values for data: min, mean, max, stdev (unbiased)""" 3950 3951 if not data: return 0,0,0,0 3952 length = len(data) 3953 if length < 1: return 0,0,0,0 3954 3955 if type(data[0]) == str: 3956 try: dd = [float(val) for val in data] 3957 except: 3958 print('** bad data for min_mean_max_stdev') 3959 return 0, 0, 0, 0 3960 else: dd = data 3961 if length == 1: return dd[0], dd[0], dd[0], 0.0 3962 3963 minval = min(dd) 3964 maxval = max(dd) 3965 meanval = loc_sum(dd)/float(length) 3966 3967 return minval, meanval, maxval, stdev_ub(dd) 3968 3969def interval_offsets(times, dur): 3970 """given a list of times and an interval duration (e.g. TR), return 3971 the offsets into respective intervals""" 3972 3973 if not times or dur <= 0: 3974 print("** interval offsets: bad dur (%s) or times: %s" % (dur, times)) 3975 return [] 3976 3977 length = len(times) 3978 if length < 1: return [] 3979 3980 fdur = float(dur) # to make sure (e.g. avoid int division) 3981 3982 try: offlist = [val % fdur for val in times] 3983 except: 3984 print("** interval offsets 2: bad dur (%s) or times: %s" % (dur, times)) 3985 return [] 3986 3987 return offlist 3988 3989def fractional_offsets(times, dur): 3990 """given a list of times and an interval duration (e.g. TR), return 3991 the fractional offsets into respective intervals 3992 3993 i.e. similar to interval offsets, but times are divided by dur""" 3994 3995 # rely on i_o for error checking 3996 olist = interval_offsets(times, dur) 3997 if len(olist) < 1 or dur <= 0: return [] 3998 3999 dur = float(dur) 4000 for ind, val in enumerate(olist): 4001 olist[ind] = val/dur 4002 4003 return olist 4004 4005def stdev_ub(data): 4006 """unbiased standard deviation (divide by len-1, not just len) 4007 stdev_ub = sqrt( (sumsq - N*mean^2)/(N-1) ) 4008 """ 4009 4010 length = len(data) 4011 if length < 2: return 0.0 4012 4013 meanval = loc_sum(data)/float(length) 4014 # compute standard deviation 4015 ssq = 0.0 4016 for val in data: ssq += val*val 4017 val = (ssq - length*meanval*meanval)/(length-1.0) 4018 4019 # watch for truncation artifact 4020 if val < 0.0 : return 0.0 4021 return math.sqrt(val) 4022 4023def stdev(data): 4024 """(biased) standard deviation (divide by len, not len-1)""" 4025 4026 length = len(data) 4027 if length < 2: return 0.0 4028 4029 meanval = loc_sum(data)/float(length) 4030 # compute standard deviation 4031 ssq = 0.0 4032 for val in data: ssq += val*val 4033 val = (ssq - length*meanval*meanval)/length 4034 4035 # watch for truncation artifact 4036 if val < 0.0 : return 0.0 4037 return math.sqrt(val) 4038 4039def variance_ub(data): 4040 """unbiased variance (divide by len-1, not just len)""" 4041 4042 length = len(data) 4043 if length < 2: return 0.0 4044 4045 meanval = loc_sum(data)/float(length) 4046 # compute standard deviation 4047 ssq = 0.0 4048 for val in data: ssq += val*val 4049 val = (ssq - length*meanval*meanval)/(length-1.0) 4050 4051 # watch for truncation artifact 4052 if val < 0.0 : return 0.0 4053 return val 4054 4055def variance(data): 4056 """(biased) variance (divide by len, not len-1)""" 4057 4058 length = len(data) 4059 if length < 2: return 0.0 4060 4061 meanval = loc_sum(data)/float(length) 4062 # compute standard deviation 4063 ssq = 0.0 4064 for val in data: ssq += val*val 4065 val = (ssq - length*meanval*meanval)/length 4066 4067 # watch for truncation artifact 4068 if val < 0.0 : return 0.0 4069 return val 4070 4071def covary(x, y): 4072 """return the raw covariance: 4073 sum[(xi - mean_x)*(yi - mean_y)] / (N-1) 4074 """ 4075 4076 ll = len(x) 4077 mx = mean(x) 4078 my = mean(y) 4079 4080 vv = loc_sum([(x[i]-mx)*(y[i]-my) for i in range(ll)]) 4081 4082 if ll > 1: return 1.0 * vv / (ll-1.0) 4083 else: return 0.0 4084 4085def r(vA, vB, unbiased=0): 4086 """return Pearson's correlation coefficient 4087 4088 for demeaned and unit vectors, r = dot product 4089 for unbiased correlation, return r*(1 + (1-r^2)/2N) 4090 4091 note: correlation_p should be faster 4092 """ 4093 la = len(vA) 4094 4095 if len(vB) != la: 4096 print('** r (correlation): vectors have different lengths') 4097 return 0.0 4098 ma = mean(vA) 4099 mb = mean(vB) 4100 dA = [v-ma for v in vA] 4101 dB = [v-mb for v in vB] 4102 sA = L2_norm(dA) # is float 4103 sB = L2_norm(dB) 4104 dA = [v/sA for v in dA] 4105 dB = [v/sB for v in dB] 4106 4107 r = dotprod(dA,dB) 4108 4109 if unbiased: return r * (1.0 + (1-r*r)/(2.0*la)) 4110 return r 4111 4112def linear_fit(vy, vx=None): 4113 """return slope and intercept for linear fit to data 4114 4115 if vx is not provided (i.e. no scatterplot fit), then return 4116 fit to straight line (i.e. apply as vx = 1..N, demeaned) 4117 4118 slope = N*sum(xy) - (sum(x)*sum(y)] 4119 --------------------------- 4120 N*sum(x*x) - (sum(x))^2 4121 4122 inter = 1/N * (sum(y) - slope*sum(x)) 4123 4124 note: we could compute slope = covary(x,y)/covary(x,x) 4125 """ 4126 4127 N = len(vy) 4128 mn = (N-1.0)/2 4129 4130 # maybe use demeaned, slope 1 line 4131 if vx == None: vx = [i-mn for i in range(N)] 4132 else: 4133 if len(vx) != N: 4134 print('** cannot fit vectors of different lengths') 4135 return 0.0, 0.0 4136 4137 sx = loc_sum(vx) 4138 sy = loc_sum(vy) 4139 ssx = dotprod(vx, vx) 4140 ssxy = dotprod(vy, vx) 4141 4142 slope = (1.0 * N * ssxy - sx * sy) / (N * ssx - sx*sx ) 4143 inter = 1.0 * (sy - slope * sx) / N 4144 4145 return slope, inter 4146 4147 4148def eta2(vA, vB): 4149 """return eta^2 (eta squared - Cohen, NeuroImage 2008 4150 4151 SUM[ (a_i - m_i)^2 + (b_i - m_i)^2 ] 4152 eta^2 = 1 - ------------------------------------ 4153 SUM[ (a_i - M )^2 + (b_i - M )^2 ] 4154 4155 where a_i and b_i are the vector elements 4156 m_i = (a_i + b_i)/2 4157 M = mean across both vectors 4158 4159 """ 4160 4161 length = len(vA) 4162 if len(vB) != length: 4163 print('** eta2: vectors have different lengths') 4164 return 0.0 4165 if length < 1: return 0.0 4166 4167 ma = mean(vA) 4168 mb = mean(vB) 4169 gm = 0.5*(ma+mb) 4170 4171 vmean = [(vA[i]+vB[i])*0.5 for i in range(length)] 4172 4173 da = [vA[i] - vmean[i] for i in range(length)] 4174 db = [vB[i] - vmean[i] for i in range(length)] 4175 num = sumsq(da) + sumsq(db) 4176 4177 da = [vA[i] - gm for i in range(length)] 4178 db = [vB[i] - gm for i in range(length)] 4179 denom = sumsq(da) + sumsq(db) 4180 4181 if num < 0.0 or denom <= 0.0 or num >= denom: 4182 print('** bad eta2: num = %s, denom = %s' % (num, denom)) 4183 return 0.0 4184 return 1.0 - float(num)/denom 4185 4186def correlation_p(vA, vB, demean=1, unbiased=0): 4187 """return the Pearson correlation between the 2 vectors 4188 (allow no demean for speed) 4189 """ 4190 4191 la = len(vA) 4192 if len(vB) != la: 4193 print('** correlation_pearson: vectors have different lengths') 4194 return 0.0 4195 4196 if la < 2: return 0.0 4197 4198 if demean: 4199 ma = mean(vA) 4200 mb = mean(vB) 4201 dA = [v-ma for v in vA] 4202 dB = [v-mb for v in vB] 4203 else: 4204 dA = vA 4205 dB = vB 4206 4207 sAB = dotprod(dA, dB) 4208 ssA = sumsq(dA) 4209 ssB = sumsq(dB) 4210 4211 if demean: del(dA); del(dB) 4212 4213 if ssA <= 0.0 or ssB <= 0.0: return 0.0 4214 else: 4215 r = sAB/math.sqrt(ssA*ssB) 4216 if unbiased: return r * (1 + (1-r*r)/(2*la)) 4217 return r 4218 4219def ttest(data0, data1=None): 4220 """just a top-level function""" 4221 4222 if data1: return ttest_2sam(data0, data1) 4223 return ttest_1sam(data0) 4224 4225def ttest_1sam(data, m0=0.0): 4226 """return (mean-m0) / (stdev_ub(data)/sqrt(N)), 4227 4228 where stdev_ub = sqrt( (sumsq - N*mean^2)/(N-1) ) 4229 4230 or faster, return: (sum-N*m0)/(sqrt(N)*stdev_ub) 4231 4232 note: move 1/N factor to denominator 4233 """ 4234 4235 # check for short length 4236 N = len(data) 4237 if N < 2: return 0.0 4238 4239 # check for division by 0 4240 sd = stdev_ub(data) 4241 if sd <= 0.0: return 0.0 4242 4243 # and return, based on any passed expected mean 4244 if m0: t = (loc_sum(data) - N*m0)/(math.sqrt(N)*sd) 4245 else: t = loc_sum(data) /(math.sqrt(N)*sd) 4246 4247 return t 4248 4249def ttest_paired(data0, data1): 4250 """easy: return 1 sample t-test of the difference""" 4251 4252 N0 = len(data0) 4253 N1 = len(data1) 4254 if N0 < 2 or N1 < 2: return 0.0 4255 if N0 != N1: 4256 print('** ttest_paired: unequal vector lengths') 4257 return 0.0 4258 4259 return ttest_1sam([data1[i] - data0[i] for i in range(N0)]) 4260 4261def ttest_2sam(data0, data1, pooled=1): 4262 """if not pooled, return ttest_2sam_unpooled(), otherwise 4263 4264 return (mean1-mean0)/sqrt(PV * (1/N0 + 1/N1)) 4265 4266 where PV (pooled_variance) = ((N0-1)*V0 + (N1-1)*V1)/(N0+N1-2) 4267 4268 note: these lists do not need to be of the same length 4269 note: the sign is as with 3dttest (second value(s) minus first) 4270 """ 4271 4272 if not pooled: return ttest_2sam_unpooled(data0, data1) 4273 4274 N0 = len(data0) 4275 N1 = len(data1) 4276 if N0 < 2 or N1 < 2: return 0.0 4277 4278 m0 = loc_sum(data0)/float(N0) 4279 v0 = variance_ub(data0) 4280 4281 m1 = loc_sum(data1)/float(N1) 4282 v1 = variance_ub(data1) 4283 4284 pv = ((N0-1)*v0 + (N1-1)*v1) / (N0+N1-2.0) 4285 if pv <= 0.0: return 0.0 4286 4287 return (m1-m0)/math.sqrt(pv * (1.0/N0 + 1.0/N1)) 4288 4289def ttest_2sam_unpooled(data0, data1): 4290 """return (mean1-mean0)/sqrt(var0/N0 + var1/N1) 4291 4292 note: these lists do not need to be of the same length 4293 note: the sign is as with 3dttest (second value(s) minus first) 4294 """ 4295 4296 N0 = len(data0) 4297 N1 = len(data1) 4298 if N0 < 2 or N1 < 2: return 0.0 4299 4300 m0 = loc_sum(data0)/float(N0) 4301 v0 = variance_ub(data0) 4302 4303 m1 = loc_sum(data1)/float(N1) 4304 v1 = variance_ub(data1) 4305 4306 if v0 <= 0.0 or v1 <= 0.0: return 0.0 4307 4308 return (m1-m0)/math.sqrt(v0/N0 + v1/N1) 4309 4310 4311def p2q(plist, do_min=1, verb=1): 4312 """convert list of p-value to a list of q-value, where 4313 q_i = minimum (for m >= i) of N * p_m / m 4314 if do min is not set, simply compute q-i = N*p_i/i 4315 4316 return q-values in increasing significance 4317 (i.e. as p goes large to small, or gets more significant) 4318 """ 4319 4320 q = plist[:] 4321 q.sort() 4322 N = len(q) 4323 4324 # work from index N down to 0 (so index using i-1) 4325 min = 1 4326 for i in range(N,0,-1): 4327 ind = i-1 4328 q[ind] = N * q[ind] / i 4329 if do_min: 4330 if q[ind] < min: min = q[ind] 4331 if min < q[ind]: q[ind] = min 4332 4333 # and flip results 4334 q.reverse() 4335 4336 return q 4337 4338def argmax(vlist, absval=0): 4339 if len(vlist) < 2: return 0 4340 if absval: vcopy = [abs(val) for val in vlist] 4341 else: vcopy = vlist 4342 4343 mval = vcopy[0] 4344 mind = 0 4345 for ind, val in enumerate(vlist): 4346 if val > mval: 4347 mval = val 4348 mind = ind 4349 4350 return mind 4351 4352def argmin(vlist, absval=0): 4353 if len(vlist) < 2: return 0 4354 if absval: vcopy = [abs(val) for val in vlist] 4355 else: vcopy = vlist 4356 4357 mval = vcopy[0] 4358 mind = 0 4359 for ind, val in enumerate(vlist): 4360 if val < mval: 4361 mval = val 4362 mind = ind 4363 4364 return mind 4365 4366def gaussian_at_hwhm_frac(frac): 4367 """gaussian_at_hwhm_frac(frac): 4368 4369 return the magnitude of a (unit, 0-centered) Gaussian curve at a 4370 fractional offset of the HWHM (half width at half maximum) = FWHM/2 4371 4372 return h(f) = 2^[-f^2] 4373 4374 HWHM is a logical input, as it is a radius, while FWHM is a diameter. 4375 4376 So gaussian_at_hwhm_frac(1) = 0.5, by definition. 4377 4378 - the return value should be in (0,1], and == 1 @ 0 4379 - the return value should be 0.5 @ 0.5 (half max @ FWHM radius) 4380 (i.e. since FWHM is a diameter, FWHM/2 is the radius) 4381 - if frac < 0, whine and return 0 (it is undefined) 4382 4383 Gaussian curves have the form: G(x) = a*e^-[ (x-b)^2 / (2*c^2) ], 4384 where a = scalar, maybe 1/(c*sqrt(2*PI)), to be unit integral 4385 b = expected value, the central point of symmetry 4386 c = standard deviation 4387 4388 A unit height, zero-centered curve has a=1, b=0: g(x)=e^-[x^2 / (2*c^2)] 4389 4390 To find (the HWHM) where g(x) = 1/2, solve: g(w) = 1/2 for w. 4391 4392 w = sqrt(c^2 * 2*ln(2)) {just use positive} 4393 4394 We want an equation for g(x), but where x is a fraction of the HWHM. 4395 Rewrite g(x) in terms of w, by solving the above for c: 4396 4397 c = w / sqrt(2 * ln2) 4398 4399 and substitute back into g(x). g(x) = e^-[x^2 * ln2 / w^2] 4400 4401 and finally write in terms of f, where f is a fraction of w. 4402 4403 Define h(f) = g(fw) = e^-[f^2*w^2 * ln2 / w^2] = e^-[f^2 * ln2] 4404 4405 Now we can cancel the 2 and ln: 4406 4407 h(f) = e^-ln[2^(f^2)] = e^ln[2^(-f^2)] = 2^(-f^2) 4408 4409 return h(f) = 2^(-f^2) 4410 """ 4411 if frac < 0: 4412 print("** gaussian_at_hwhm_frac: illegal frac < 0 of %s", frac) 4413 return 0 4414 4415 return 2.0 ** -(frac*frac) 4416 4417def gaussian_at_fwhm(x, fwhm): 4418 """gaussian_at_fwhm(x, fwhm): 4419 4420 return the magnitude of unit, zero-centered Gaussian curve at the given x 4421 4422 The Gaussian curve is defined by the given FWHM value, e.g. 4423 4424 g(x) = e^-[x^2 * ln2 / FWHM^2] 4425 4426 This actually returns gaussian_at_hwhm_frac(x/(fwhm/2)), or of (2x/fwhm). 4427 """ 4428 if fwhm <= 0: 4429 print("** gaussian_at_fwhm: illegal fwhm <= 0 of %s", fwhm) 4430 return 0 4431 4432 return gaussian_at_hwhm_frac(2.0*x/fwhm) 4433 4434# ---------------------------------------------------------------------- 4435# random list routines: shuffle, merge, swap, extreme checking 4436# ---------------------------------------------------------------------- 4437 4438def swap_2(vlist, i1, i2): 4439 if i1 != i2: 4440 val = vlist[i2] 4441 vlist[i2] = vlist[i1] 4442 vlist[i1] = val 4443 4444def shuffle(vlist, start=0, end=-1): 4445 """randomize the order of list elements, where each permutation is 4446 equally likely 4447 4448 - akin to RSFgen, but do it with equal probabilities 4449 (search for swap in [index,N-1], not in [0,N-1]) 4450 - random.shuffle() cannot produce all possibilities, don't use it 4451 - start and end are indices to work with 4452 """ 4453 4454 # if we need random elsewhere, maybe do it globally 4455 import random 4456 4457 vlen = len(vlist) 4458 4459 # check bounds 4460 if start < 0 or start >= vlen: return 4461 if end >= 0 and end <= start: return 4462 4463 # so start is valid and end is either < 0 or big enough 4464 if end < 0 or end >= vlen: end = vlen-1 4465 4466 nvals = end-start+1 4467 4468 # for each index, swap with random other towards end 4469 for index in range(nvals-1): 4470 rind = int((nvals-index)*random.random()) 4471 swap_2(vlist, start+index, start+index+rind) 4472 continue 4473 4474 # return list reference, though usually ignored 4475 return vlist 4476 4477def shuffle_blocks(vlist, bsize=-1): 4478 """like shuffle, but in blocks of given size""" 4479 4480 vlen = len(vlist) 4481 4482 if bsize < 0 or bsize >= vlen: 4483 shuffle(vlist) 4484 return 4485 4486 if bsize < 2: return 4487 4488 nblocks = vlen // bsize 4489 nrem = vlen % bsize 4490 4491 boff = 0 4492 for ind in range(nblocks): 4493 shuffle(vlist, boff, boff+bsize-1) 4494 boff += bsize 4495 shuffle(vlist, boff, boff+nrem-1) 4496 4497 return vlist 4498 4499def random_merge(list1, list2): 4500 """randomly merge 2 lists (so each list stays in order) 4501 4502 shuffle a list of 0s and 1s and then fill from lists 4503 """ 4504 4505 # if we need random elsewhere, maybe do it globally 4506 import random 4507 4508 mlist = [0 for i in range(len(list1))] + [1 for i in range(len(list2))] 4509 shuffle(mlist) 4510 4511 i1, i2 = 0, 0 4512 for ind in range(len(mlist)): 4513 if mlist[ind] == 0: 4514 mlist[ind] = list1[i1] 4515 i1 += 1 4516 else: 4517 mlist[ind] = list2[i2] 4518 i2 += 1 4519 4520 return mlist 4521 4522def show_sum_pswr(nT, nR): 4523 cp = 0.0 4524 prev = 0 4525 for r in range(nR+1): 4526 # already float, but be clear 4527 p = float(prob_start_with_R(nT,nR,r)) 4528 cp += p 4529 # print 'prob at %3d = %g (cum %g)' % (r, p, cp) 4530 if prev == 0: prev = p 4531 print(p, p/prev) 4532 prev = p 4533 print('cum result is %g' % cp) 4534 4535 4536def prob_start_with_R(nA, nB, nS): 4537 """return the probability of starting nS (out of nB) class B elements 4538 should equal: choose(nB, nS)*nS! * nA *(nB+nA-nS-1)! / (nA+nB)! 4539 or: factorial(nB, init=nB-nS+1) * nA / fact(nA+nB, init=nA+nB-nS) 4540 4541 or: choose(nB,nS)/choose(nA+nB,nS) * nA/(nA+nB-nS) 4542 4543 """ 4544 return 1.0 * nA * factorial(nB, init=nB-nS+1) \ 4545 / factorial(nA+nB, init=nA+nB-nS) 4546 4547def choose(n,m): 4548 """return n choose m = n! / (m! * (n-m)!)""" 4549 # integral division (or use floats, to better handle overflow) 4550 return factorial(n,init=n-m+1) / factorial(m) 4551 4552def factorial(n, init=1): 4553 prod = 1 4554 val = init 4555 while val <= n: 4556 prod *= val 4557 val += 1 4558 return prod 4559 4560def swap2(data): 4561 """swap data elements in pairs""" 4562 4563 size = 2 4564 nsets = len(data)//size 4565 if nsets <= 0: return 4566 4567 for ind in range(nsets): 4568 off = ind*size 4569 v = data[off] # swap d[0] and d[1] 4570 data[off] = data[off+1] 4571 data[off+1] = v 4572 4573def swap4(data): 4574 """swap data elements in groups of 4""" 4575 4576 size = 4 4577 nsets = len(data)//size 4578 if nsets <= 0: return 4579 4580 for ind in range(nsets): 4581 off = ind*size 4582 v = data[off] # swap d[0] and d[3] 4583 data[off] = data[off+3] 4584 data[off+3] = v 4585 v = data[off+1] # swap d[1] and d[2] 4586 data[off+1] = data[off+2] 4587 data[off+2] = v 4588 4589def vec_extremes(vec, minv, maxv, inclusive=0): 4590 """return a integer array where values outside bounds are 1, else 0 4591 4592 if inclusive, values will also be set if they equal a bound 4593 4594 return error code, new list 4595 success: 0, list 4596 error : 1, None""" 4597 4598 if not vec: return 1, None 4599 4600 if minv > maxv: 4601 print('** extremes: minv > maxv (', minv, maxv, ')') 4602 return 1, None 4603 4604 if inclusive: 4605 elist = [1*(vec[t]>=maxv or vec[t]<=minv) for t in range(len(vec))] 4606 else: 4607 elist = [1*(vec[t]> maxv or vec[t]< minv) for t in range(len(vec))] 4608 4609 return 0, elist 4610 4611def vec_moderates(vec, minv, maxv, inclusive=1): 4612 """return a integer array where values inside bounds are 1, else 0 4613 4614 if inclusive, values will also be set if they equal a bound 4615 4616 return error code, new list 4617 success: 0, list 4618 error : 1, None""" 4619 4620 if not vec: return 1, None 4621 4622 if minv > maxv: 4623 print('** moderates: minv > maxv (', minv, maxv, ')') 4624 return 1, None 4625 4626 if inclusive: 4627 elist = [1*(vec[t]>=minv and vec[t]<=maxv) for t in range(len(vec))] 4628 else: 4629 elist = [1*(vec[t]> minv and vec[t]< maxv) for t in range(len(vec))] 4630 4631 return 0, elist 4632 4633def vec_range_limit(vec, minv, maxv): 4634 """restrict the values to [minv, maxv] 4635 4636 This function modifies the past vector. 4637 4638 return 0 on success, 1 on error""" 4639 4640 if not vec: return 0 4641 4642 if minv > maxv: 4643 print('** range_limit: minv > maxv (', minv, maxv, ')') 4644 return 1 4645 4646 for ind in range(len(vec)): 4647 if vec[ind] < minv: vec[ind] = minv 4648 elif vec[ind] > maxv: vec[ind] = maxv 4649 4650 return 0 4651 4652# for now, make 2 vectors and return their correlation 4653def test_polort_const(ntrs, nruns, verb=1): 4654 """compute the correlation between baseline regressors of length ntrs*nruns 4655 - make vectors of 11...10000...0 and 00...011...100..0 that are as the 4656 constant polort terms of the first 2 runs 4657 - return their correlation 4658 4659 - note that the correlation is easily provable as -1/(N-1) for N runs 4660 """ 4661 4662 if ntrs <= 0 or nruns <= 2: return -1 # flag 4663 4664 # lazy way to make vectors 4665 v0 = [1] * ntrs + [0] * ntrs + [0] * (ntrs * (nruns-2)) 4666 v1 = [0] * ntrs + [1] * ntrs + [0] * (ntrs * (nruns-2)) 4667 4668 if verb > 1: 4669 print('++ test_polort_const, vectors are:\n' \ 4670 ' v0 : %s \n' \ 4671 ' v1 : %s' % (v0, v1)) 4672 4673 return correlation_p(v0, v1) 4674 4675# for now, make 2 vectors and return their correlation 4676def test_tent_vecs(val, freq, length): 4677 a = [] 4678 b = [] 4679 for i in range(length): 4680 if (i%freq) == 0: 4681 a.append(val) 4682 b.append(1-val) 4683 elif ((i-1)%freq) == 0: 4684 a.append(0.0) 4685 b.append(val) 4686 else: 4687 a.append(0.0) 4688 b.append(0.0) 4689 4690 return correlation_p(a,b) 4691 4692 4693# ----------------------------------------------------------------------- 4694# [PT: June 8, 2020] for matching str entries in list (for the FATCAT 4695# -> MVM and other group analysis programs) 4696## [PT: June 8, 2020] updated to allow cases where only a subset of 4697## either A or B is matched; but will still give errors if something 4698## matches more than 1 item in the other list. 4699 4700def match_listA_str_in_listB_str(A, B): 4701 """Input: two lists (A and B), each of whose elements are strings. A 4702 and B don't have to have the same length. 4703 4704 See if each string in A is contained is contained within one (and 4705 only one) string in list B. If yes, return: 4706 1 4707 the dictionary of matches, A->B 4708 the dictionary of matches, B->A 4709 elif only a subset of A is contained in B or vice versa, return: 4710 0 4711 the dictionary of matches, A->B 4712 the dictionary of matches, B->A 4713 otherwise, error+exit. 4714 4715 The primary/first usage of this program is for the following case: 4716 matching subject IDs from a CSV file (the list of which form A) 4717 with path names for *.grid|*.netcc files (the list of which form 4718 B), for combining all that information. We want to find 1 subject 4719 ID from the CSV file with 1 (and only 1) matrix file, for each 4720 subject ID. 4721 4722 NB: in the above case, the users are alternatively able to enter 4723 the list for matching subj ID with matrix file, if the above 4724 name-matching wouldn't work. 4725 4726 """ 4727 4728 if type(A) != list or type(B) != list : 4729 BASE.EP("Both inputs A and B must be lists, not {} and {}, " 4730 "respectively".format(type(A), type(B))) 4731 4732 na = len(A) 4733 nb = len(B) 4734 4735 if not(na and nb) : 4736 BASE.EP("One of the sets is empty: len(A) = {}; len(B) = {}" 4737 "".format(na, nb)) 4738 4739 # check that each ele is a str 4740 ta = [type(x)!=str for x in A] 4741 tb = [type(x)!=str for x in B] 4742 if max(ta) or max(tb) : 4743 BASE.EP("All elements of A and B must be str, but that is not true\n" 4744 "for a least one list: for A, it's {}; for B it's {}" 4745 "".format(not(max(ta)), not(max(tb)))) 4746 4747 checklist_a = [0]*na 4748 checklist_b = [0]*nb 4749 matchlist_a = [-1]*na 4750 matchlist_b = [-1]*nb 4751 4752 for ii in range(nb): 4753 for jj in range(na): 4754 if B[ii].__contains__(A[jj]) : 4755 checklist_a[jj]+= 1 4756 checklist_b[ii]+= 1 4757 matchlist_a[jj] = ii 4758 matchlist_b[ii] = jj 4759 # *should* be able to break here, but will check 4760 # *all*, to verify there are no problems/ambiguities 4761 4762 CHECK_GOOD = True # determines if we 'return' anything 4763 FULL_MATCH = 1 # flags type of matching when returning 4764 4765 # --------- now check the outcomes 4766 if min(checklist_a) == 1 and max(checklist_a) == 1 : 4767 # all matches found, all singletons 4768 BASE.IP("Found single matches for each element in A") 4769 else: 4770 FULL_MATCH = 0 4771 if min(checklist_a) == 0 : 4772 # all found matches are singletons, but there are gaps; 4773 # not a fatal error 4774 BASE.WP("Some elements of A are unmatched:") 4775 for ii in range(na): 4776 if not(checklist_a[ii]) : 4777 print("\t unmatched: {}".format(A[ii])) 4778 if max(checklist_a) > 1 : 4779 # ambiguities in matching --> badness 4780 CHECK_GOOD = False 4781 BASE.WP("Some elements of A are overmatched:") 4782 for ii in range(na): 4783 if checklist_a[ii] > 1 : 4784 print("\t overmatched: {}".format(A[ii])) 4785 if not(max(checklist_a)) : 4786 CHECK_GOOD = False 4787 BASE.WP("No elements in A matched??:") 4788 4789 if min(checklist_b) == 1 and max(checklist_b) == 1 : 4790 # all matches found, all singletons 4791 BASE.IP("Found single matches for each element in B") 4792 else: 4793 FULL_MATCH = 0 4794 if min(checklist_b) == 0 : 4795 # all found matches are singletons, but there are gaps; 4796 # not a fatal error 4797 BASE.WP("Some elements of B are unmatched:") 4798 for ii in range(nb): 4799 if not(checklist_b[ii]) : 4800 print("\t unmatched: {}".format(B[ii])) 4801 if max(checklist_b) > 1 : 4802 # ambiguities in matching --> badness 4803 CHECK_GOOD = False 4804 BASE.WP("Some elements of B are overmatched:") 4805 for ii in range(nb): 4806 if checklist_b[ii] > 1 : 4807 print("\t overmatched: {}".format(B[ii])) 4808 if not(max(checklist_b)) : 4809 CHECK_GOOD = False 4810 BASE.WP("No elements in B matched??:") 4811 4812 if not(CHECK_GOOD) : 4813 BASE.EP('Exiting') 4814 4815 # if we made it here, things are good 4816 da = {} 4817 for ii in range(na): 4818 if checklist_a[ii] : 4819 da[ii] = matchlist_a[ii] 4820 db = {} 4821 for ii in range(nb): 4822 if checklist_b[ii] : 4823 db[ii] = matchlist_b[ii] 4824 4825 na_keys = len(da.keys()) 4826 nb_keys = len(db.keys()) 4827 4828 # final check on consistency 4829 if na_keys != nb_keys : 4830 BASE.EP("Mismatch in number of keys in output lists?\n" 4831 "There are {} for A and {} for B" 4832 "".format(na_keys, nb_keys)) 4833 4834 # e.g., 4835 #for key in da: 4836 # print(A[key] + ' <---> ' + B[da[key]]) 4837 #for key in db: 4838 # print(A[db[key]] + ' <---> ' + B[key]) 4839 4840 return FULL_MATCH, da, db 4841 4842def invert_dict(D): 4843 """Take a dictionary D of key-value pairs and return a dictionary 4844 Dinv, which the keys of D as values and the (matched) values of D 4845 as keys. 4846 4847 """ 4848 4849 if not(D) : return {} 4850 4851 Dinv = {} 4852 for k, v in D.items(): 4853 Dinv[v] = k 4854 4855 return Dinv 4856 4857 4858# ---------------------------------------------------------------------- 4859# [PT: Jan 13, 2020] Pieces for dealing with files of AFNI seed points 4860 4861class afni_seeds: 4862 4863 def __init__(self, file_line=[]): 4864 '''Make an object for seed-based correlation. 4865 4866 Use the optional 'file_line' input to read a (non-split) line 4867 from a file easily. 4868 4869 ''' 4870 4871 self.xyz = [] 4872 self.space = '' 4873 self.roi_label = '' 4874 self.netw = None 4875 4876 # these for future/other use; cannot be changed at the moment 4877 self.rad = None 4878 self.shape = 'sphere' 4879 4880 # use this to read a (non-split) line from a file easily 4881 if file_line : 4882 self.set_all_info_from_file(file_line) 4883 4884 def set_space(self, x): 4885 self.space = x 4886 4887 def set_xyz(self, xlist): 4888 '''Input a list of strings or floats to be a list of float coords. 4889 ''' 4890 if len(xlist) != 3: 4891 print("** ERROR: need 3 items input, not {}:\n" 4892 " '{}'".format(len(xlist), xlist)) 4893 sys.exit(3) 4894 self.xyz = [float(x) for x in xlist] 4895 4896 def set_roi_label(self, x): 4897 self.roi_label = x 4898 4899 def set_netw(self, x): 4900 self.netw = x 4901 4902 def set_rad(self, x): 4903 self.rad = float(x) 4904 4905 def set_shape(self, x): 4906 self.shape = x 4907 4908 def set_all_info_from_file(self, line): 4909 '''Read in an unsplit line from a file, and save all the information 4910 appropriately. Column order matters! 4911 4912 ''' 4913 y = line.split() 4914 Ny = len(y) 4915 if Ny < 5 : 4916 print("** ERROR: too few columns ({}) in this line:\n" 4917 " {}\n" 4918 " e.g., did you forget something?".format(Ny, line)) 4919 sys.exit(2) 4920 4921 elif Ny > 6 : 4922 print("** ERROR: too many columns ({}) in this line:\n" 4923 " {}\n" 4924 " e.g., did you put spaces in names?".format(Ny, line)) 4925 sys.exit(2) 4926 4927 xyz = [ float(y[j]) for j in range(3) ] 4928 4929 self.set_xyz(xyz) 4930 self.set_space(y[3]) 4931 self.set_roi_label(y[4]) 4932 4933 if len(y) == 6 : 4934 self.set_netw(y[5]) 4935 4936# [PT: Jan 13, 2020] for APQC (and maybe other purposes later), read 4937# in a simple text file that contains both numbers and strings. 4938# See text file: afni_seeds_per_space.txt 4939def read_afni_seed_file(fname, only_from_space=None): 4940 '''Read in a file of columns of info, and return list of seed objs. 4941 4942 Input 4943 ----- 4944 fname : (req) text file of seed point information. 4945 only_from_space : (opt) can choose to only attach seeds from a given space 4946 (which is entered here); otherwise, all seeds are output. 4947 4948 Output: list of seed objects. 4949 4950 ''' 4951 4952 # expected/allowed cols ('netw' is opt) 4953 cols = ['x', 'y', 'z', 'space', 'roi_label', 'netw'] 4954 Ncols = len(cols) 4955 4956 # input data, ignore blank lines 4957 x = read_text_file( fname, noblank=1 ) 4958 Nx = len(x) 4959 4960 dat = [] 4961 4962 for i in range(Nx): 4963 row = x[i] 4964 if row[0] != '#': 4965 seed_obj = afni_seeds(file_line=row) 4966 if not(only_from_space) : 4967 dat.append( seed_obj ) 4968 else: 4969 if seed_obj.space == only_from_space : 4970 dat.append( seed_obj ) 4971 4972 return dat 4973 4974# ---------------------------------------------------------------------- 4975 4976if __name__ == '__main__': 4977 print('afni_util.py: not intended as a main program') 4978 print(' (consider afni_python_wrapper.py)') 4979 sys.exit(1) 4980 4981