1#!/usr/local/bin/python3.8 2 3# python3 status: compatible 4 5import sys, os 6import copy, glob 7from afnipy import afni_util as UTIL 8from afnipy import lib_vars_object as VO 9 10g_mema_tests = [None, 'paired', 'unpaired'] 11g_ttpp_tests = ['-AminusB', '-BminusA'] 12 13# atomic (type within nested list) and simple types for VarsObject 14g_valid_atomic_types = [int, float, str, list] 15g_simple_types = [int, float, str] 16g_subject_sort_key = None 17 18def comment_section_string(comment, length=70, cchar='-'): 19 """return a string of the form: 20 # -------------- some comment -------------- 21 where the total length is given 22 """ 23 clen = len(comment) 24 25 ndash = (length - clen - 4) // 2 # one '#' and 3 spaces 26 27 # if no space for multiple dashes, don't use any 28 if ndash < 2: return '# %s' % comment 29 30 dstr = cchar * ndash # make one ------------- string 31 32 return '# %s %s %s\n' % (dstr, comment, dstr) 33 34def make_message_list_string(mlist, title): 35 if len(mlist) == 0: return '' 36 mesg = '' 37 for ind, mm in enumerate(mlist): 38 if ind == 0: mesg += mm 39 else: mesg += ('\n' + mm) 40 return mesg 41 42 43# =========================================================================== 44# begin Subject stuff (class should be rewritten to use VarsObject types) 45# =========================================================================== 46 47def subj_compare_key(subj): 48 """subj_compare() will no longer be valid in python3, as sort() will 49 lose the cmp= attribute. Generate a key that would be equivalent, 50 and just use reverse() after the fact, if desired. 51 52 The key is now stored in g_subject_sort_key. 53 54 Return a [key, sid] pair (or sid, if key==None). 55 """ 56 if g_subject_sort_key == None: 57 return subj.sid 58 59 if g_subject_sort_key in subj.atrs: 60 val = subj.atrs[g_subject_sort_key] 61 else: 62 val = None 63 64 return [val, subj.sid] 65 66def subj_compare(subj0, subj1): 67 """compare 2 Subject objects: used for sorting a list of subjects 68 if _sort_key is set, compare by key, then by sid 69 else, compare by sid 70 """ 71 72 cval = 0 73 key = subj0._sort_key 74 if key != None: 75 if key in subj0.atrs: v0 = subj0.atrs[key] 76 else: v0 = None 77 if key in subj1.atrs: v1 = subj1.atrs[key] 78 else: v1 = None 79 cval = cmp(v0, v1) 80 81 # remove subj0._order, as compare must be replaced by key 82 if cval != 0: return cval 83 return cmp(subj0.sid, subj1.sid) 84 85def set_var_str_from_def(obj, name, vlist, vobj, defs, 86 verb=1, csort=1, spec=None): 87 """try to set name = value based on vlist 88 (just set as string) 89 if name is not known by the defaults, return failure 90 91 if spec: update_vars_from_special(csort) 92 93 This function will generally be called via a user interface library 94 or in the user interface itself. Since we are setting variables as 95 string 96 97 return 1 on change, 0 on unchanged, -1 on error 98 """ 99 100 if not defs.valid(name): 101 print('** invalid %s variable: %s' % (obj, name)) 102 return -1 103 104 dtype = type(defs.val(name)) 105 if dtype not in g_valid_atomic_types: 106 print('** SVSFD: unknown %s variable type for %s' % (obj, name)) 107 return -1 108 109 # if simple type but have list, fail 110 if dtype != list and len(vlist) > 1: 111 print("** SVSFD: simple variable '%s' %s\n" \ 112 " but have list value: %s" % (name, dtype, vlist)) 113 return -1 114 115 # ---------------------------------------- 116 # try to apply the value (list) 117 118 val = None 119 120 # update val, 121 # check that the it can be properly converted (if simple type) 122 if defs.has_simple_type(name): 123 val = vlist[0] 124 try: vv = dtype(val) 125 except: 126 print('** SVSFD %s.%s, cannot convert value %s to %s' \ 127 % (obj, name, val, dtype)) 128 return -1 129 elif dtype == list: val = vlist 130 else: 131 print('** SVSFD: invalid type %s for %s'%(dtype,name)) 132 return -1 133 134 # actually set the value 135 rv = vobj.set_var(name, val) 136 if verb > 1: 137 if rv: print('++ %s: updating %s to %s %s' % (obj, name, val, type(val))) 138 else: print('++ %s: no update for %s to %s' % (obj, name, val)) 139 140 # if no update, we're outta here 141 if rv == 0: return rv 142 143 # ---------------------------------------------------------------------- 144 # handle some special cases, such as indices and labels, which might 145 # come with file name lists 146 147 # this function must be passed, since it will vary per library 148 149 if spec != None: spec(name, vobj, check_sort=csort) 150 151 return rv 152 153def goto_proc_dir(dname): 154 """ go to processing directory, returning return_dir 155 - if proc_dir does not exist, create it 156 - cd 157 - return ret_dir 158 """ 159 if UTIL.is_trivial_dir(dname): return '' # nowhere to go 160 161 retdir = os.getcwd() # so need return directory 162 163 # if the directory does not yet exist, create it 164 if not os.path.isdir(dname): 165 try: os.makedirs(dname) 166 except: 167 print('** failed makedirs(%s)' % dname) 168 return '' 169 170 # now try to go there 171 try: os.chdir(dname) 172 except: 173 print('** failed to go to process dir, %s' % dname) 174 return '' 175 176 return retdir # only returned on success 177 178def ret_from_proc_dir(rname): 179 """if retdir is set, cd to retdir 180 (should be called 'atomically' with goto_proc_dir) 181 return '', to use to clear previous return dir""" 182 183 # if retdir is useless, bail 184 if rname == None or rname == '' or rname == '.': return '' 185 186 try: os.chdir(rname) 187 except: 188 print('** failed to return to %s from process dir' % rname) 189 190 return '' 191 192def proc_dir_file_exists(dname, fname): 193 if UTIL.is_trivial_dir(dname): pathname = fname 194 else: pathname = '%s/%s' % (dname, fname) 195 196 return os.path.isfile(pathname) 197 198def get_def_tool_path(prog_name, top_dir='tool_results', prefix='tool', 199 keep_if_missing='', digits=3): 200 """return something of the form top_dir/prefix.0001.prog_name 201 202 if top_dir exists: 203 - look for anything of the form prefix.* 204 - if keep_is_missing and the last entry does not contain it, 205 then return that last entry 206 ** goal is to only create new directory when process has happened 207 - find the lowest index that is not used 208 - return top_dir/prefix.NEW_INDEX.tname 209 210 else: return top_dir/prefix.001.tname 211 212 e.g. tool_results/tool.004.align_test 213 e.g. group_results/test.004.3dttest++ 214 """ 215 216 tname = prog_name # yeah, shorter, but less descriptive ... 217 tdir = top_dir 218 index = 1 # default index 219 220 # generate form for name (e.g. tr/t.%03d.t); insert index later 221 form = '%s/%s.%%0%dd.%s' % (tdir, prefix, digits, tname) 222 sform = '%s/%s.%%0%dd' % (tdir, prefix, digits) # short form 223 224 # if tdir does not yet exist, we can start with the default 225 if not os.path.isdir(tdir): 226 return form % index 227 228 # see what is under tdir, and go with default if nothing is found 229 glist = glob.glob('%s/%s.*.*' % (tdir, prefix)) 230 if len(glist) == 0: return form % index 231 232 # get trailing directory name 233 tlist = [name.split('/')[-1] for name in glist] 234 235 # abuse '.': make a list of integers from field 1 when split over '.' 236 try: ilist = [int(name.split('.')[1]) for name in tlist] 237 except: 238 print('** found non-int VAL in %s/%s.VAL.*' % (tdir, prefix)) 239 return form % 999 240 241 ilist.sort() 242 nvals = len(ilist) 243 244 # if the keep_if_missing file does NOT exist underneath the last dir, 245 # return the directory for that index 246 if keep_if_missing: 247 fdir = sform % ilist[-1] 248 if not UTIL.glob_form_has_match('%s.*/%s' % (fdir, keep_if_missing)): 249 return form % ilist[-1] # longer form includes tool name 250 251 # quick check, if ilist[n-1] <= n, just increment 252 # (or should we forget this, since non-unique values break logic?) 253 # (no, tool name may change without 'keep' file/dir, so values may repeat) 254 if ilist[-1] <= nvals: return form % (ilist[-1]+1) 255 256 # finally! now find first value > index+1 (else, just use next) 257 258 for ind, ival in enumerate(ilist): 259 if ival > ind+1: break 260 261 return form % (ind+1) 262 263class Subject(object): 264 """a simple subject object holding an ID, dataset name, and an 265 attribute dictionary""" 266 267 _sort_key = None # attribute key used in compare() 268 _order = 1 # 1 for normal, -1 for reverse 269 270 def __init__(self, sid='', dset='', atrs=None): 271 272 self.sid = sid # subject ID (string) 273 self.dset = dset # dset name (string: existing filename) 274 self.atrs = None # attributes (VarsObject instance) 275 276 self.ddir = '.' # split dset name into directory and file 277 self.dfile = '' 278 self.maxlinelen = 0 # if set, maximum line len for subj in command 279 # rcr - todo 280 281 dir, file = os.path.split(dset) # and update them 282 if dir: self.ddir = dir 283 self.dfile = file 284 285 # init to empty, and merge if something is passed 286 self.atrs = VO.VarsObject('subject %s' % sid) 287 if atrs != None: self.atrs.merge(atrs) 288 289 def show(self): 290 natr = self.atrs.count()-1 # do not include name 291 print("Subject %s, dset %s, natr = %d" % (self.sid, self.dset, natr)) 292 print(" ddir = %s\n dfile = %s\n" % (self.ddir, self.dfile)) 293 if natr > 0: 294 self.atrs.show(' attributes: ') 295 296class SubjectList(object): 297 """list of Subject elements, with attributes and command writing functions 298 299 - can pass list of SIDs, dset names, attributes (VarsObject instances) 300 - if any lists are past, the lengths must be equal 301 302 - given a list of datasets (and no SIDs), try to extract SIDs 303 by parsing into a regular expression 304 305 - sort comparison should be by _sort_key (else sid), then by sid 306 """ 307 308 def __init__(self, name='subject list', sid_l=None, dset_l=None, atr_l=None, 309 verb=1): 310 311 self.name = name # in case there are multiple lists in use 312 self.subjects = [] # list of Subject instances 313 self.atrl = [] # list of subject attributes (names only) 314 self.disp_atrs = [] # attributes to display, in order 315 self.common_dir = '' # common parent dir to subject dsets 316 self.common_dname = '' # variable to apply for common_dir ($ddir) 317 self.verb = verb # verbose level 318 self.status = 0 # non-zero is bad 319 320 if sid_l == None and dset_l == None and atr_l == None: return 321 322 # set the length, and check for consistency 323 llen = -1 324 errs = 0 325 for ilist in [sid_l, dset_l, atr_l]: 326 if ilist == None: continue 327 if llen < 0: llen = len(ilist) # set, or 328 elif llen != len(ilist): errs += 1 # ... test 329 330 if errs > 0: 331 print('** SubjectList init requires equal lengths (or None)') 332 print(' sid_l = %s' % sid_l) 333 print(' dset_l = %s' % dset_l) 334 print(' atr_l = %s' % atr_l) 335 self.status = 1 336 return 337 338 if llen == 0: return # empty lists? 339 340 # okay, fill the subjects and atrs fields 341 sid = '' 342 dname = '' 343 atr = None 344 for ind in range(llen): 345 if sid_l != None: sid = sid_l[ind] 346 if dset_l != None: dset = dset_l[ind] 347 if atr_l != None: atr = atr_l[ind] 348 self.add(Subject(sid=sid, dset=dset, atrs=atr)) 349 350 def copy(self, sid_l=[], atr='', atrval=None): 351 """make a full copy of the SubjectList element 352 353 start with a complete copy 354 if len(sid_l)>0, remove subjects not in list 355 if atr is not '', remove subjects for which atrs.atr!=atrval 356 """ 357 # start with everything and then delete 358 olen = len(newSL.subjects) 359 newSL = copy.deepcopy(self) 360 if len(sid_l) > 0: 361 skeep = [] 362 snuke = [] 363 for subj in self.subjects: 364 if subj.sid in sid_l: skeep.append(subj) 365 else : snuke.append(subj) 366 if self.verb > 2: 367 print('++ SL copy: nuking %d subjs not in sid list' % len(snuke)) 368 for subj in snuke: del(subj) 369 newSL.subjects = skeep 370 if atr != '': 371 skeep = [] 372 snuke = [] 373 for subj in self.subjects: 374 if subj.atrs.val(atr) == atrval: skeep.append(subj) 375 else: snuke.append(subj) 376 if self.verb>2: print('++ SL copy: nuking %d subjs with atr[%s] != %s'\ 377 % (atr, atrval)) 378 for subj in snuke: del(subj) 379 newSL.subjects = skeep 380 if self.verb > 1: print('++ SL copy: keeping %d of %d subjects' \ 381 % (len(newSL.subjects), olen)) 382 return newSL 383 384 def show(self, mesg=''): 385 if mesg: mstr = " (%s)" % mesg 386 else: mstr = '' 387 print("SubjectList: %s%s" % (self.name, mstr)) 388 print(" nsubj = %d, natrs = %d, ndisp_atrs = %d" % \ 389 (len(self.subjects), len(self.atrl), len(self.disp_atrs))) 390 print(" atrl = %s" % self.atrl) 391 print(" disp_atrs = %s" % self.disp_atrs) 392 393 if len(self.subjects) == 0: return 394 395 print(" subject sort key: %s" % self.subjects[0]._sort_key) 396 for subj in self.subjects: 397 subj.show() 398 399 def add(self, subj): 400 """add the subject to the list and update the atrl list""" 401 402 for atr in subj.atrs.attributes(): 403 if not atr in self.atrl: self.atrl.append(atr) 404 405 self.atrl.sort() 406 407 self.subjects.append(subj) 408 409 def set_common_data_dir(self, cname='data_dir'): 410 """return the directory common to all subject ddir names""" 411 cdir = UTIL.common_dir([s.dset for s in self.subjects]) 412 if UTIL.is_trivial_dir(cdir) or (len(cdir) < len(cname)): 413 self.common_dir = '' 414 self.common_dname = '' 415 else: 416 self.common_dir = cdir 417 self.common_dname = cname 418 if self.verb > 1: 419 print('++ setting common dir, %s = %s' % (cname, cdir)) 420 421 def set_ids_from_dsets(self, prefix='', suffix='', hpad=0, tpad=0, dpre=0): 422 """use the varying part of the dataset names for subject IDs 423 424 If hpad > 0 or tpad > 0, expand into the head or tail of the dsets. 425 If prefix or suffix is passed, apply them. 426 427 return 0 on success, 1 on error 428 """ 429 430 if hpad < 0 or tpad < 0: 431 print('** set_ids_from_dsets: will not apply negative padding') 432 return 1 433 434 # try filenames without paths, first 435 dlist = [s.dset.split('/')[-1] for s in self.subjects] 436 if UTIL.vals_are_constant(dlist): 437 print('** constant dataset names (%s)' % dlist[0]) 438 print(' trying directories...') 439 dlist = [s.dset for s in self.subjects] 440 441 slist = UTIL.list_minus_glob_form(dlist, hpad, tpad, keep_dent_pre=dpre) 442 443 # in the case of diretories, check for success 444 # (maybe we can try to skip past them, that might be okay) 445 for index in range(len(slist)): 446 if '/' in slist[index]: 447 posn = slist[index].rfind('/') 448 slist[index] = slist[index][posn+1:] 449 if len(slist[index]) < 1: 450 print('** failed to extract subject IDs from directory list') 451 print(' (directories do not vary at single level)') 452 return 1 453 454 if len(slist) != len(self.subjects): 455 print('** failed to set SIDs from dset names\n' \ 456 ' dsets = %s\n' \ 457 ' slist = %s' % (dlist, slist)) 458 return 1 459 460 if not UTIL.vals_are_unique(slist): 461 print('** cannot set IDs from dsets, labels not unique: %s' % slist) 462 print('-- labels come from dsets: %s' % dlist) 463 return 1 464 465 for ind, subj in enumerate(self.subjects): 466 subj.sid = '%s%s%s' % (prefix, slist[ind], suffix) 467 468 return 0 469 470 def restrict_ids_to_dsets(self, valid_ids=[], require=1): 471 """restrict subject IDs to those in valid_ids list 472 require all valid_ids to exist, or fail 473 474 return 0 on success 475 """ 476 # bail if either list is empty 477 if len(self.subjects) == 0: return 0 478 if len(valid_ids) == 0: return 0 479 480 # check that valid_ids are unique 481 if not UTIL.vals_are_unique(valid_ids): 482 print('** restrict_ids: ids are not unique') 483 return 1 484 485 # check that all valid_ids exist, and generate new subject list 486 all_ids = [subj.sid for subj in self.subjects] 487 488 new_subjs = [] 489 missing = 0 490 missed_id = '' # example of missing ID 491 for sid in valid_ids: 492 if sid in all_ids: 493 old_index = all_ids.index(sid) 494 new_subjs.append(self.subjects[old_index]) 495 else: 496 if self.verb > 1: 497 print("** restrict_ids: cannot restrict to missing ID '%s'"%sid) 498 missed_id = sid 499 missing += 1 500 if missing: 501 print("** restrict_ids: missing %d of %d IDs" \ 502 % (missing,len(valid_ids))) 503 print(" IDs look like: %s" % ' '.join(all_ids[:3])) 504 print(" missing IDs look like: %s" % missed_id) 505 if require: 506 return 1 507 else: 508 print("-- restrict_ids: allowing %d missing IDs..." % missing) 509 510 # apply restricted list 511 self.subjects = new_subjs 512 513 return 0 514 515 def remove_ids_from_dsets(self, remove_ids=[], require=1): 516 """restrict subject IDs to those not in remove_ids list 517 if require: require all remove_ids to exist, or fail 518 519 return 0 on success 520 """ 521 # bail if either list is empty 522 if len(self.subjects) == 0: return 0 523 if len(remove_ids) == 0: return 0 524 525 # check that remove_ids are unique 526 if not UTIL.vals_are_unique(remove_ids): 527 print('** remove_ids: ids are not unique') 528 return 1 529 530 # check that all remove_ids exist, and fail if not 531 all_ids = [subj.sid for subj in self.subjects] 532 missing = 0 533 for sid in remove_ids: 534 if sid not in all_ids: 535 if self.verb > 1: 536 print("** remove_ids: cannot remove missing ID '%s'"%sid) 537 missed_id = sid 538 missing += 1 539 if missing and (require or self.verb > 1): 540 print("** remove_ids: missing %d of %d IDs" \ 541 % (missing,len(remove_ids))) 542 print(" IDs look like: %s" % ' '.join(all_ids[:3])) 543 print(" missing IDs look like: %s" % missed_id) 544 # if required, this is fatal 545 if require: 546 return 1 547 else: 548 print("-- remove_ids: allowing %d missing IDs..." % missing) 549 550 # generate a new subject list 551 new_subjs = [] 552 for sindex, sid in enumerate(all_ids): 553 if sid not in remove_ids: 554 new_subjs.append(self.subjects[sindex]) 555 556 # apply remove list 557 self.subjects = new_subjs 558 559 return 0 560 561 def sort(self, key=None, order=1): 562 if len(self.subjects) == 0: return 563 # sort() has no cmp keyword in python3, use key method 564 # Subject._sort_key = key # None or otherwise 565 # Subject._order = order # 1 for small first, -1 for reverse 566 # self.subjects.sort(cmp=subj_compare) 567 568 g_subject_sort_key = key 569 self.subjects.sort(key=subj_compare_key) 570 if order < 0: self.subjects.reverse() 571 572 def make_anova2_command(self, bsubs=None, prefix=None, options=None, verb=1): 573 """create a basic 3dANOVA2 -type 3 command 574 575 ** bsubs should be lists of strings, even if integral sub-bricks 576 (they are applied as sub-brick selectors) 577 578 attach options after subject lists 579 580 bsubs - beta sub-bricks (1 list of sub-brick selectors) 581 prefix - prefix for command output 582 options - other options added to the command 583 verb - verbose level 584 585 return None on failure, command on success 586 """ 587 588 if prefix == '' or prefix == None: prefix = 'anova2_result' 589 if verb > 1: print('-- make_anova2_command: have prefix %s' % prefix) 590 591 if bsubs == None: 592 print('** missing sub-brick selection list') 593 return None 594 if len(bsubs) < 2: 595 print('** anova2_command: need at least 2 sub-bricks (have %d)' \ 596 % len(bsubs)) 597 return None 598 599 indent = 9 # minimum indent: spaces to following -set option 600 601 cmd = '#!/bin/tcsh\n\n' 602 603 # maybe we will use directory variables 604 self.set_common_data_dir() 605 if not UTIL.is_trivial_dir(self.common_dir): 606 self.common_dname = 'data' 607 cmd += '# apply any data directories with variables\n' \ 608 'set %s = %s\n' % (self.common_dname, self.common_dir) 609 cmd += '\n' 610 611 cmd += '# note: factor A is condition, B is subject\n\n' 612 613 # command and first set of subject files 614 cmd += '3dANOVA2 -type 3 \\\n' \ 615 '%s' % self.make_anova2_set_list(bsubs, indent) 616 617 if len(options) > 0: cmd += '%*s%s \\\n' % (indent,'', ' '.join(options)) 618 else: # add some basic option 619 opt = '-amean 1 amean1' 620 cmd += '%*s%s \\\n' % (indent, '', opt) 621 print('++ no contrast options given, adding simple: %s' % opt) 622 623 if prefix.find('/') >= 0: pp = prefix 624 else: pp = './%s' % prefix 625 cmd += '%*s-bucket %s\n' % (indent, '', pp) 626 627 cmd += '\n' 628 629 return cmd 630 631 def make_anova3_command(self, bsubs=None, prefix=None, subjlists=None, 632 options=None, factors=[], verb=1): 633 """create a basic 3dANOVA3 -type 5 command 634 635 ** other types may be added later... 636 637 ** bsubs should be lists of strings, even if integral sub-bricks 638 (they are applied as sub-brick selectors) 639 640 attach options after subject lists 641 642 bsubs - beta sub-bricks (1 list of sub-brick selectors) 643 prefix - prefix for command output 644 subjlists - len > 1 for type 5 645 options - other options added to the command 646 atype - 3dANOVA3 -type (should be 4 or 5) 647 factors - if type 4, #factors of each type (f0*f1 = len(b)) 648 verb - verbose level 649 650 Note: for type 5: factor A is group, B is condition, C is subject 651 652 return None on failure, command on success 653 """ 654 655 if prefix == '' or prefix == None: prefix = 'anova3_result' 656 if verb > 1: print('-- make_anova2_command: have prefix %s' % prefix) 657 658 if bsubs == None: 659 print('** missing sub-brick selection list') 660 return None 661 if len(bsubs) < 2: 662 print('** anova3_command: need at least 2 sub-bricks (have %d)' \ 663 % len(bsubs)) 664 return None 665 666 ncond = len(factors) 667 ngroups = len(subjlists) 668 669 atype = 0 670 if ngroups > 1: atype = 5 671 elif ncond == 2: atype = 4 672 673 if atype == 4: 674 if ngroups != 1: 675 print('** anova3_cmd: -type 4 requires only 1 dset group') 676 return None 677 if ncond != 2: 678 print('** anova3_cmd: -type 4 requires 2 factor lengths') 679 print(' (product should be length -subs_betas)') 680 return None 681 if factors[0]*factors[1] != len(bsubs): 682 print('** anova3_cmd: -type 4 factor mismatch') 683 print(' (%d x %d design requires %d betas, have %d' \ 684 % (factors[0], factors[1], factors[0]*factors[1], len(bsubs))) 685 return None 686 elif atype == 5: 687 if ngroups < 2: 688 print('** anova3_cmd: -type 5 requires >= 2 subject lists') 689 return None 690 if ncond > 1: 691 print('** anova3_cmd: -type 5 should not have sets of factors') 692 return None 693 else: 694 print('** anova3_cmd: cannot detect -type 4 or -type 5, seek -help!') 695 return None 696 697 indent = 4 # indent after main command 698 699 # maybe we will use directory variables 700 cmd = '#!/bin/tcsh\n\n' 701 found = 0 702 slen0 = len(subjlists[0].subjects) 703 subjlists[0].set_common_data_dir() 704 cd0 = subjlists[0].common_dir 705 706 for ilist, slist in enumerate(subjlists): 707 slist.set_common_data_dir() 708 if not UTIL.is_trivial_dir(slist.common_dir): 709 if not found: # first time found 710 cmd += '# apply any data directories with variables\n' 711 found = 1 712 if ilist > 0 and slist.common_dir == cd0: 713 slist.common_dname = 'data1' 714 else: 715 slist.common_dname = 'data%d' % (ilist+1) 716 cmd += 'set %s = %s\n' % (slist.common_dname, slist.common_dir) 717 if found: cmd += '\n' 718 719 if atype == 4: 720 cmd += '# note: factor A is cond 1, B is cond 2, C is subject\n\n' 721 sfunc = self.make_anova3_t4_set_list 722 else: 723 cmd += '# note: factor A is group, B is condition, C is subject\n\n' 724 sfunc = self.make_anova3_t5_set_list 725 726 # command and first set of subject files 727 cstr = sfunc(bsubs, subjlists, factors, indent) 728 if cstr == None: return None 729 cmd += '3dANOVA3 -type %d \\\n' '%s' % (atype, cstr) 730 731 if len(options) > 0: cmd += '%*s%s \\\n' % (indent,'', ' '.join(options)) 732 else: # add some basic options 733 opt = '-amean 1 amean1 -bmean 1 bmean1' 734 cmd += '%*s-amean 1 amean1 \\\n' \ 735 '%*s-bmean 1 bmean1 \\\n' % (indent, '', indent, '') 736 print('++ no contrast options given, adding simple: %s' % opt) 737 738 if prefix.find('/') >= 0: pp = prefix 739 else: pp = './%s' % prefix 740 cmd += '%*s-bucket %s\n' % (indent, '', pp) 741 742 cmd += '\n' 743 744 return cmd 745 746 def make_generic_command(self, command, bsubs=None, subjlist2=None, 747 prefix=None, options=''): 748 """create a generic command 749 750 This basically allows one to create a generic command that takes 751 AFNI-style inputs. 752 753 Note: this is almost identical to make_mema_command, except for use 754 of tsubs. 755 756 attach options before subject lists 757 bsubs can be of length 1 even with 2 set_labs, in that case: 758 - length 1: must have subjlist2 set, and apply to it 759 - length 2: no subjlist2, use with subjlist 1 760 761 return None on failure, command on success 762 """ 763 764 if not command: return 765 766 # do we have a second slist? 767 s2 = subjlist2 # so much typing... 768 769 # make sure we have sub-brick selectors for possibly both lists 770 if bsubs == None: bsubs = [None] 771 if s2 and len(bsubs) == 1: bsubs.append(bsubs[0]) 772 773 # ready for work, maybe note status 774 if self.verb > 1: 775 print('-- make_generic_command: %s -prefix %s' % (command, prefix)) 776 print(' s2 = %d, bsubs = %s'%(s2!=None,bsubs)) 777 778 # initialize directories and variables 779 rv, cmd = self.set_data_dirs(subjlist2=s2) 780 if rv: return 781 782 # append actual command, prefix and options 783 indent = len(command)+1 784 785 # having no prefix is valid, skip that part 786 if prefix: cmd += '%s -prefix %s \\\n' % (command, prefix) 787 else: cmd += '%s \\\n' % (command) 788 789 if len(options) > 0: 790 cmd += '%*s%s \\\n' % (indent, '', ' '.join(options)) 791 792 # if s2: cmd += '%*s# dataset list 1 \\\n' % (indent,'') 793 cmd += self.make_generic_set_list(bsubs[0], indent) 794 if s2: 795 cmd += '\\\n' + s2.make_generic_set_list(bsubs[1], indent) 796 elif len(bsubs) > 1: 797 cmd += '\\\n' + self.make_generic_set_list(bsubs[1], indent) 798 799 # strip trailing backslash (must dupe memory) 800 if cmd[-2:] == '\\\n': cmd = cmd[0:-2] 801 802 cmd += '\n\n' 803 804 return cmd 805 806 def set_data_dirs(self, subjlist2=None): 807 """Given 1 or 2 file lists, set directories and return initial 808 script to apply them with variables. 809 810 return status and script text 811 812 if there are 0 or 1 common dirs use them 813 if there are 2 and they are different, adjust var names 814 """ 815 s2 = subjlist2 # so much typing... 816 self.set_common_data_dir() 817 uses2dir = 0 818 diffdirs = 0 819 if s2 != None: 820 s2.set_common_data_dir() 821 uses2dir = not UTIL.is_trivial_dir(s2.common_dir) 822 if not UTIL.is_trivial_dir(s2.common_dir) and \ 823 not UTIL.is_trivial_dir(self.common_dir) and \ 824 self.common_dir != s2.common_dir: 825 # differentiate the directory variable names 826 self.common_dname = 'data1' 827 s2.common_dname = 'data2' 828 diffdirs = 1 829 830 cmd = '' 831 if not UTIL.is_trivial_dir(self.common_dir) or uses2dir: 832 cmd += '# apply any data directories with variables\n' 833 834 if not UTIL.is_trivial_dir(self.common_dir): 835 cmd += 'set %s = %s\n' % (self.common_dname, self.common_dir) 836 837 if diffdirs: 838 cmd += 'set %s = %s\n' % (s2.common_dname, s2.common_dir) 839 840 if cmd: cmd += '\n' 841 842 return 0, cmd 843 844 845 def make_ttestpp_command(self, set_labs=None, bsubs=None, subjlist2=None, 846 prefix=None, comp_dir=None, options=None, verb=1): 847 """create a basic 3dttest++ command 848 849 Note: this is almost identical to make_mema_command, except for use 850 of tsubs. 851 852 ** labs, bsubs should be lists of strings, even if they are integral 853 sub-bricks 854 855 if set_labs=None, use defaults depending on # of subject lists 856 else, set_labs must have 1 or 2 elements, for 1 or 2 sets of subjects 857 bsubs can be of length 1 even with 2 set_labs, in that case: 858 - length 1: must have subjlist2 set, and apply to it 859 - length 2: no subjlist2, use with subjlist 1 860 attach options after subject lists 861 862 set_labs - set labels (None, 1 or 2 labels) 863 bsubs - beta sub-bricks (length 1 or 2, matching labels) 864 subjlist2 - second subject list for 2-sample test (want if the 865 datasets differ across sets) 866 prefix - prefix for 3dtest++ output 867 comp_dir - comparison direction, either -AminusB or -BminusA 868 (if 2 sets) 869 options - other options added to the 3dtest++ command 870 verb - verbose level 871 872 return None on failure, command on success 873 """ 874 875 if prefix == '' or prefix == None: prefix = 'ttest++_result' 876 if verb > 1: print('-- make_ttest++_command: have prefix %s' % prefix) 877 s2 = subjlist2 # sooooo much typing... 878 879 if set_labs == None: 880 if s2 == None: set_labs = ['setA'] 881 else: set_labs = ['setA', 'setB'] 882 if verb > 2: print('-- tt++_cmd: adding default set labels') 883 if bsubs == None: bsubs, tsubs = ['0'], ['1'] 884 885 indent = 3 # minimum indent: spaces to following -set option 886 887 # want any '-AminusB' option at top of command 888 if len(set_labs) > 1: copt = '%*s%s \\\n' % (indent, '', comp_dir) 889 else: copt = '' 890 891 # maybe we will use directory variables 892 self.set_common_data_dir() 893 if not UTIL.is_trivial_dir(self.common_dir): 894 if s2 != None: 895 s2.set_common_data_dir() 896 if UTIL.is_trivial_dir(s2.common_dir): 897 # then do not use either 898 self.common_dir = '' 899 self.common_dname = '' 900 # else, use both (if same, default variable is okay) 901 elif self.common_dir != s2.common_dir: 902 # different, so update the names to be different 903 self.common_dname = 'data1' 904 s2.common_dname = 'data2' 905 906 cmd = '' 907 if not UTIL.is_trivial_dir(self.common_dir): 908 cmd += '# apply any data directories with variables\n' \ 909 'set %s = %s\n' % (self.common_dname, self.common_dir) 910 if s2 != None: 911 if not UTIL.is_trivial_dir(s2.common_dir) \ 912 and s2.common_dir != self.common_dir: 913 cmd += 'set %s = %s\n' % (s2.common_dname, s2.common_dir) 914 cmd += '\n' 915 916 # command and first set of subject files 917 cmd += '3dttest++ \\\n' \ 918 '%*s-prefix %s \\\n' \ 919 '%s' \ 920 '%*s-setA %s \\\n%s' \ 921 % (indent, ' ', prefix, copt, 922 indent, ' ', set_labs[0], 923 self.make_ttpp_set_list(bsubs[0], indent+3)) 924 925 # maybe add second set of subject files 926 if len(set_labs) > 1: 927 if verb > 2: print('-- tt++_cmd: have labels for second subject set') 928 929 # separate tests for bad test types 930 if comp_dir == None: 931 print('** make_tt++_cmd: missing test type, should be in %s' \ 932 % g_ttpp_tests) 933 return # failure 934 if comp_dir not in g_ttpp_tests: 935 print('** make_tt++_cmd: comp_dir (%s) must be in the list %s' \ 936 % (comp_dir, g_ttpp_tests)) 937 return # failure 938 939 # note subject list and sub-brick labels 940 if s2 != None: 941 S = s2 942 if verb > 2: print('-- second subject list was passed') 943 else: 944 S = self 945 if verb > 2: print('-- no second subject list, using same list') 946 if len(bsubs) > 1: b = bsubs[1] 947 else: 948 if S != s2: 949 print('** make_tt++_cmd: same subject list in comparison') 950 b = bsubs[0] 951 cmd += '%*s-setB %s \\\n%s' % \ 952 (indent, ' ', set_labs[1], S.make_ttpp_set_list(b, indent+3)) 953 954 if len(options) > 0: cmd += '%*s%s' % (indent, '', ' '.join(options)) 955 956 # strip trailing backslash (must dupe memory) 957 if cmd[-2:] == '\\\n': cmd = cmd[0:-2] 958 959 cmd += '\n\n' 960 961 return cmd 962 963 def make_mema_command(self, set_labs=None, bsubs=None, tsubs=None, 964 subjlist2=None, prefix=None, ttype=None, options=None, 965 verb=1): 966 """create a basic 3dMEMA command 967 968 ** labs, bsubs, tsubs should be lists of strings, even if they are 969 integral sub-bricks 970 971 if set_labs=None, use defaults depending on # of subject lists 972 else, set_labs must have 1 or 2 elements, for 1 or 2 sets of subjects 973 bsubs and tsubs can be of length 1 even with 2 set_labs, in that case: 974 - length 1: must have subjlist2 set, and apply to it 975 - length 2: no subjlist2, use with subjlist 1 976 attach options after subject lists 977 978 set_labs - set labels (None, 1 or 2 labels) 979 bsubs - beta sub-bricks (length 1 or 2, matching labels) 980 tsubs - t-stat sub-bricks (as with bsubs) 981 subjlist2 - second subject list for 2-sample test (want if the 982 datasets differ across sets) 983 prefix - prefix for 3dMEMA output 984 ttype - used for a 2-sample test, to distinguish between 985 paired and unpaired (results is using either 986 -conditions (paired) or -group (un-)) 987 ** ttype == paired is no longer valid 988 options - other options added to the 3dMEMA command 989 verb - verbose level 990 991 return None on failure, command on success 992 """ 993 994 if prefix == '' or prefix == None: prefix = 'mema_result' 995 if verb > 1: print('++ make_mema_cmd: have prefix %s' % prefix) 996 s2 = subjlist2 # sooooo much typing... 997 998 if set_labs == None: 999 if s2 == None: set_labs = ['setA'] 1000 else: set_labs = ['setA', 'setB'] 1001 if verb > 2: print('++ mema_cmd: adding default set labels') 1002 if bsubs == None: bsubs, tsubs = ['0'], ['1'] 1003 1004 # maybe we will use directory variables 1005 self.set_common_data_dir() 1006 if not UTIL.is_trivial_dir(self.common_dir): 1007 if s2 != None: 1008 s2.set_common_data_dir() 1009 if UTIL.is_trivial_dir(s2.common_dir): 1010 # then do not use either 1011 self.common_dir = '' 1012 self.common_dname = '' 1013 # else, use both (if same, default variable is okay) 1014 elif self.common_dir != s2.common_dir: 1015 # different, so update the names to be different 1016 self.common_dname = 'data1' 1017 s2.common_dname = 'data2' 1018 1019 cmd = '' 1020 if not UTIL.is_trivial_dir(self.common_dir): 1021 cmd += '# apply any data directories with variables\n' \ 1022 'set %s = %s\n' % (self.common_dname, self.common_dir) 1023 if s2 != None: 1024 if not UTIL.is_trivial_dir(s2.common_dir) \ 1025 and s2.common_dir != self.common_dir: 1026 cmd += 'set %s = %s\n' % (s2.common_dname, s2.common_dir) 1027 cmd += '\n' 1028 1029 # command and first set of subject files 1030 cmd += '3dMEMA -prefix %s \\\n' \ 1031 ' -set %s \\\n%s' % \ 1032 (prefix,set_labs[0],self.make_mema_set_list(bsubs[0],tsubs[0],10)) 1033 1034 # maybe add second set of subject files 1035 if len(set_labs) > 1: 1036 if verb > 2: print('-- mema_cmd: have labels for second subject set') 1037 # note subject list and sub-brick labels 1038 if s2 != None: 1039 S = s2 1040 if verb > 2: print('-- second subject list was passed') 1041 else: 1042 S = self 1043 if verb > 2: print('-- no second subject list, using same list') 1044 if len(bsubs) > 1: b, t = bsubs[1], tsubs[1] 1045 else: 1046 if S != s2: 1047 print('** make_mema_cmd: same subject list in comparison') 1048 b, t = bsubs[0], tsubs[0] 1049 cmd += '%7s-set %s \\\n%s' % \ 1050 (' ', set_labs[1], S.make_mema_set_list(b, t, 10)) 1051 1052 # either 2-sample or paired 1053 if ttype not in g_mema_tests: 1054 print("** invalid 3dMEMA test %s, not in %s" % (ttype,g_mema_tests)) 1055 return None 1056 if ttype == 'paired': 1057 print('** 3dMEMA -type paired: no longer valid\n' \ 1058 ' (input contrast and t-stat from original regression)') 1059 return None 1060 else: opt = '-groups' 1061 cmd += '%7s%s %s %s \\\n' % ('', opt, set_labs[0], set_labs[1]) 1062 1063 if len(options) > 0: cmd += '%7s%s' % ('', ' '.join(options)) 1064 1065 # strip trailing backslash (must dupe memory) 1066 if cmd[-2:] == '\\\n': cmd = cmd[0:-2] 1067 1068 cmd += '\n\n' 1069 1070 return cmd 1071 1072 def make_mema_set_list(self, bsub, tsub, indent=0): 1073 """return a multi-line string of the form: 1074 SID1 "dset1[bsub]" 1075 "dset1[tsub]" 1076 SID2 "dset2[bsub]" 1077 "dset2[tsub]" 1078 ... 1079 indent is the initial indentation 1080 """ 1081 # note the max subject ID length 1082 ml = 0 1083 for subj in self.subjects: 1084 if len(subj.sid) > ml: ml = len(subj.sid) 1085 1086 if not UTIL.is_trivial_dir(self.common_dir) and self.common_dname: 1087 sdir = self.common_dname 1088 else: sdir = '' 1089 1090 sstr = '' 1091 for subj in self.subjects: 1092 if sdir: 1093 # see if the dataset is in a directory underneath 1094 cdir = UTIL.child_dir_name(self.common_dir, subj.ddir) 1095 if UTIL.is_trivial_dir(cdir): cstr = '' 1096 else: cstr = '%s/' % cdir 1097 dset = '$%s/%s%s' % (sdir, cstr, subj.dfile) 1098 else: dset = subj.dset 1099 sstr += '%*s%s "%s[%s]" \\\n%*s "%s[%s]" \\\n' % \ 1100 (indent, '', subj.sid, dset, bsub, 1101 indent+ml, '', dset, tsub) 1102 return sstr 1103 1104 def make_generic_set_list(self, bsub, indent=0): 1105 """return a multi-line string of the form: 1106 "dset1[bsub]" 1107 "dset2[bsub]" 1108 ... 1109 indent is per-line indentation 1110 """ 1111 if not UTIL.is_trivial_dir(self.common_dir) and self.common_dname: 1112 sdir = self.common_dname 1113 else: sdir = '' 1114 1115 sstr = '' 1116 for subj in self.subjects: 1117 if sdir: 1118 # see if the dataset is in a directory underneath 1119 cdir = UTIL.child_dir_name(self.common_dir, subj.ddir) 1120 if UTIL.is_trivial_dir(cdir): cstr = '' 1121 else: cstr = '%s/' % cdir 1122 dset = '$%s/%s%s' % (sdir, cstr, subj.dfile) 1123 else: dset = subj.dset 1124 if bsub == None: 1125 sstr += '%*s%s \\\n' % (indent, '', dset) 1126 else: # use bsub 1127 sstr += '%*s"%s[%s]" \\\n' % (indent, '', dset, bsub) 1128 1129 return sstr 1130 1131 def make_ttpp_set_list(self, bsub, indent=0): 1132 """return a multi-line string of the form: 1133 SID1 "dset1[bsub]" 1134 SID2 "dset2[bsub]" 1135 ... 1136 indent is the initial indentation 1137 """ 1138 # note the max subject ID length 1139 ml = 0 1140 for subj in self.subjects: 1141 if len(subj.sid) > ml: ml = len(subj.sid) 1142 1143 if not UTIL.is_trivial_dir(self.common_dir) and self.common_dname: 1144 sdir = self.common_dname 1145 else: sdir = '' 1146 1147 sstr = '' 1148 for subj in self.subjects: 1149 if sdir: 1150 # see if the dataset is in a directory underneath 1151 cdir = UTIL.child_dir_name(self.common_dir, subj.ddir) 1152 if UTIL.is_trivial_dir(cdir): cstr = '' 1153 else: cstr = '%s/' % cdir 1154 dset = '$%s/%s%s' % (sdir, cstr, subj.dfile) 1155 else: dset = subj.dset 1156 sstr += '%*s%s "%s[%s]" \\\n' % (indent, '', subj.sid, dset, bsub) 1157 1158 return sstr 1159 1160 def make_anova2_set_list(self, bsub, indent=0): 1161 """return a multi-line string of the form: 1162 -alevels #bsub 1163 -blevels #subj 1164 -dset ALEVEL BLEVEL "dset#A[bsub#B]" 1165 ... 1166 indent is the initial indentation 1167 """ 1168 sdir = self.common_dname 1169 sstr = '%*s-alevels %d \\\n' \ 1170 '%*s-blevels %d \\\n' \ 1171 % (indent,'', len(bsub), indent, '', len(self.subjects)) 1172 1173 for isubj, subj in enumerate(self.subjects): 1174 if sdir: 1175 # see if the dataset is in a directory underneath 1176 cdir = UTIL.child_dir_name(self.common_dir, subj.ddir) 1177 if UTIL.is_trivial_dir(cdir): cstr = '' 1178 else: cstr = '%s/' % cdir 1179 dset = '$%s/%s%s' % (sdir, cstr, subj.dfile) 1180 else: dset = subj.dset 1181 for ibeta, beta in enumerate(bsub): 1182 sstr += '%*s-dset %2d %2d "%s[%s]" \\\n' \ 1183 % (indent, '', ibeta+1, isubj+1, dset, beta) 1184 1185 return sstr 1186 1187 def make_anova3_t4_set_list(self, bsub, subjlists, factors, indent=0): 1188 """return a multi-line string of the form: 1189 -alevels #alevels 1190 -blevels #blevels 1191 -clevels #subj 1192 -dset ALEVEL BLEVEL SUBJ "dset#A[bsub#B]" 1193 ... 1194 - factors should be of length 2 1195 - indent is the initial indentation 1196 - as in type 5, have A change slower than B, but subj be slowest 1197 (so factor order per subject matches command line) 1198 """ 1199 1200 errs = 0 1201 1202 if len(subjlists) != 1: 1203 print('** MAt4SL: bad subject list count = %d' % len(subjlists)) 1204 return None 1205 1206 nA = factors[0] 1207 nB = factors[1] 1208 slist = subjlists[0] 1209 if nA*nB != len(bsub): 1210 print('** MAt4SL: bad factor count: %d, %d, %d' % (nA, nB, len(bsub))) 1211 return None 1212 1213 sstr = '' 1214 sstr += '%*s-alevels %d \\\n' % (indent, '', nA) 1215 sstr += '%*s-blevels %d \\\n' % (indent, '', nB) 1216 sstr += '%*s-clevels %d \\\n' % (indent, '', len(slist.subjects)) 1217 1218 sdir = slist.common_dname 1219 for isubj, subj in enumerate(slist.subjects): 1220 if sdir: 1221 # see if the dataset is in a directory underneath 1222 cdir = UTIL.child_dir_name(slist.common_dir, subj.ddir) 1223 if UTIL.is_trivial_dir(cdir): cstr = '' 1224 else: cstr = '%s/' % cdir 1225 dset = '$%s/%s%s' % (sdir, cstr, subj.dfile) 1226 else: dset = subj.dset 1227 1228 for iA in range(nA): 1229 for iB in range(nB): 1230 sstr += '%*s-dset %2d %2d %2d "%s[%s]" \\\n' \ 1231 % (indent, '', iA+1, iB+1, isubj+1, dset, 1232 bsub[iA*nB+iB]) 1233 1234 if errs: return None 1235 1236 return sstr 1237 1238 def make_anova3_t5_set_list(self, bsub, subjlists, factors=0, indent=0): 1239 """return a multi-line string of the form: 1240 -alevels #subjlists 1241 -blevels #bsub 1242 -clevels #subj 1243 -dset GROUP BLEVEL SUBJ "dset#A[bsub#B]" 1244 ... 1245 factors is ignored, and exists only to match type4 function 1246 indent is the initial indentation 1247 """ 1248 sstr = '' 1249 sstr += '%*s-alevels %d \\\n' % (indent, '', len(subjlists)) 1250 sstr += '%*s-blevels %d \\\n' % (indent, '', len(bsub)) 1251 sstr += '%*s-clevels %d \\\n' % (indent, '', len(subjlists[0].subjects)) 1252 1253 slen0 = len(subjlists[0].subjects) 1254 errs = 0 1255 for ilist, slist in enumerate(subjlists): 1256 if len(slist.subjects) != slen0: 1257 print('** subject list %d length differs from SL 1 (%d != %d)\n'\ 1258 % (ilist+1, len(slist.subjects), slen0)) 1259 errs += 1 1260 sdir = slist.common_dname 1261 1262 for isubj, subj in enumerate(slist.subjects): 1263 if sdir: 1264 # see if the dataset is in a directory underneath 1265 cdir = UTIL.child_dir_name(slist.common_dir, subj.ddir) 1266 if UTIL.is_trivial_dir(cdir): cstr = '' 1267 else: cstr = '%s/' % cdir 1268 dset = '$%s/%s%s' % (sdir, cstr, subj.dfile) 1269 else: dset = subj.dset 1270 for ibeta, beta in enumerate(bsub): 1271 sstr += '%*s-dset %2d %2d %2d "%s[%s]" \\\n' \ 1272 % (indent, '', ilist+1, ibeta+1, isubj+1, dset, beta) 1273 1274 if errs: return None 1275 1276 return sstr 1277 1278if __name__ == '__main__': 1279 print('** this is not a main program') 1280 1281