1#!/usr/local/bin/python3.8 2 3# python3 status: started 4 5# currently, this explicitly does _not_ depend on scipy or numpy 6 7from __future__ import print_function 8 9 10import os, sys 11from afnipy import module_test_lib 12g_testlibs = ['math', 'copy'] 13if module_test_lib.num_import_failures(g_testlibs): sys.exit(1) 14 15 16# import libraries 17import math 18import copy 19 20from afnipy import afni_util as UTIL 21from afnipy import afni_base as BASE 22from afnipy import lib_textdata as TD 23 24MTYPE_NONE = 0 # no modulation (these work as a bit mask) 25MTYPE_AMP = 1 # amplitude modulation 26MTYPE_DUR = 2 # duration modulation 27 28# global variables 29g_xmat_basis_labels = ['Name', 'Option', 'Formula', 'Columns'] 30g_xmat_stim_types = [ 'times', 'AM', 'AM1', 'AM2', 'IM' ] 31 32class Afni1D: 33 def __init__(self, filename="", from_mat=0, matrix=None, verb=1): 34 """1D matrix class (2-D array) 35 36 matrix is stored transposed from 1D file, each row as a time series 37 38 init from filename or matrix, 39 filename : 1D file to read in as matrix 40 from_mat : use matrix for initialization ([]) 41 matrix : 2-D matrix to use if from_mat is set 42 name : some text to associate with the datase 43 verb : verbose level of operations 44 45 mat : data (transposed, as vectors of data over time) 46 fname : filename 47 nvec, nt : dimensions 48 labels : simple text labels 49 groups : list: -1 (baseline), 0 (motion), else reg index 50 tr : in seconds 51 """ 52 53 # main variables 54 self.mat = None # actual data (2-D array [[]]) 55 self.name = "NoName" # more personal and touchy-feely... 56 self.fname = filename # name of data file 57 self.aname = None # afni_name, if parsing is useful 58 59 # .xmat.1D variables 60 self.nvec = 0 # one vector per column in file 61 self.nt = 0 # length of each column (# rows) 62 self.tr = 1.0 63 self.nrowfull = 0 # length without censoring 64 self.nruns = 1 65 self.run_len = [0] # len(run_len) is number of runs 66 self.run_len_nc= [0] # run lengths without censoring 67 self.nroi = 1 68 self.GLapplied = 0 # was goodlist applied (1=yes, 2=zero-padded) 69 70 # misc variables (from attributes) 71 self.command = '' # from CommandLine 72 self.header = [] # array of extra header (comment) lines 73 self.csimobj = None # ClustSim object 74 self.csimstat = -1 # -1, 0, 1: undef, not csim, is csim 75 76 # list variables (from attributes) 77 self.havelabs = 0 # did we find any labels 78 self.labels = [] # label per vector (from ColumnLabels) 79 self.groups = [] # integral column groups (from ColumnGroups) 80 self.goodlist = [] # good time points (from GoodList) 81 self.runstart = [] # start indices (from RunStart) 82 self.basisinfo= [] # Basis* fields as dictionaries 83 # Name, Option, Formula, Columns 84 85 self.verb = verb 86 self.ready = 0 # matrix is ready 87 88 # computed variables 89 self.cormat = None # correlation mat (normed xtx) 90 self.cosmat = None # cosine mat (scaled xtx) 91 self.cormat_ready = 0 # correlation mat is set 92 93 self.VO = None # VarsObject, if set 94 95 # initialize... 96 if self.fname: 97 if self.init_from_general_name(self.fname): return None 98 elif from_mat: 99 if self.init_from_matrix(matrix): return None 100 101 self.update_group_info() 102 103 if self.verb > 2: self.show() 104 105 def reduce_by_tlist(self, tlist): 106 """reduce by time list, similiar to afni's {} selector 107 this affects run_len and runs, so try to fix those 108 109 return 0 on success""" 110 111 if not self.ready: 112 print('** append: Afni1D is not ready') 113 return 1 114 if not UTIL.is_valid_int_list(tlist, 0, self.nt-1, whine=1): return 1 115 if len(tlist) < 1: return 1 116 117 self.mat = [[row[t] for t in tlist] for row in self.mat] 118 119 if self.nrowfull == self.nt: self.nrowfull = len(tlist) 120 else: self.nrowfull = 0 # cannot tell 121 self.nt = len(tlist) 122 123 if self.goodlist: self.goodlist = [self.goodlist[t] for t in tlist] 124 125 self.update_group_info() # since nruns or run_len may have changed 126 127 return 0 128 129 def reduce_by_run_list(self, rlist): 130 """extract the list of (1-based) runs from the dataset 131 return 0 on success""" 132 133 if not self.ready: 134 print('** reduce_by_run_list: Afni1D is not ready') 135 return 1 136 137 # are all of the runs valid? 138 errs = 0 139 for run in rlist: 140 if run < 1: 141 print('** reduce_by_run_list: min run index is 1') 142 errs += 1 143 if run > self.nruns: 144 print('** reduce_by_run_list: %d exceeds max run index %d' \ 145 % (run, self.nruns)) 146 errs += 1 147 if errs: return 1 148 149 if len(self.runstart) != self.nruns: 150 print('** run_list: missing %d run positions (-set_run_lengths?)' \ 151 % self.nruns) 152 errs += 1 153 if len(self.run_len) != self.nruns: 154 print('** run_list: missing %d run lengths (-set_run_lengths?)' \ 155 % self.nruns) 156 errs += 1 157 if errs: return 1 158 159 # make useful runstart list to begin with, which will be trashed anyway 160 runstart = [] 161 ntotal = 0 162 for rlen in self.run_len: 163 runstart.append(ntotal) 164 ntotal += rlen 165 self.runstart = runstart 166 167 # finally ready? 168 # make tlist, clear groups, reduce_by_tlist 169 # and then update nruns, run_len, runstart 170 self.groups = [] 171 tlist = [] 172 run_len = [] # make new lists as we go 173 runstart = [] # make new lists as we go 174 ntotal = 0 # for runstart list 175 for run in rlist: 176 r0 = run - 1 177 nt = self.run_len[r0] 178 tlist.extend([self.runstart[r0]+i for i in range(nt)]) 179 run_len.append(nt) 180 runstart.append(ntotal) 181 ntotal += nt 182 183 if self.reduce_by_tlist(tlist): return 1 184 185 self.nruns = len(rlist) 186 self.run_len = run_len 187 self.runstart = runstart 188 189 return 0 190 191 def reduce_by_vec_list(self, vlist): 192 """reduce the dataset according to the vector list 193 return 0 on success""" 194 195 if not self.ready: 196 print('** reduce_by_vec_list: Afni1D is not ready') 197 return 1 198 if len(vlist) == 0: 199 # allow creation of empty matrices, but possibly whine 200 if self.verb > 1: 201 print("** afni1D: reducing to empty matrix") 202 elif not UTIL.is_valid_int_list(vlist, 0, self.nvec-1, whine=1): 203 return 1 204 205 self.mat = [self.mat[i] for i in vlist] 206 207 self.nvec = len(vlist) 208 self.nroi = self.nvec # don't know, so assume all 209 210 if self.nvec == 0: self.nt = 0 211 212 # update lists 213 if self.labels: self.labels = [self.labels[i] for i in vlist] 214 if self.groups: self.groups = [self.groups[i] for i in vlist] 215 216 if self.basisinfo: 217 # consider propagating this later 218 del(self.basisinfo) 219 self.basisinfo = [] 220 221 return 0 222 223 def reduce_by_group_list(self, glist): 224 """reduce the dataset according to the list of groups 225 return 0 on success""" 226 227 if not self.ready: 228 print('** reduce_by_group_list: Afni1D is not ready') 229 return 1 230 231 glen = len(glist) 232 if glen < 1: return 0 233 234 if len(self.groups) < 1: 235 print('** no groups to select in %s' % self.name) 236 return 1 237 238 # do not all repeats in group list 239 groups = UTIL.get_unique_sublist(self.groups) 240 241 gnew = glist[:] 242 243 # convert to an integral group list 244 for ind in range(glen-1,-1,-1): 245 if gnew[ind] == 'POS': 246 gnew[ind:ind+1] = [g for g in groups if g > 0] 247 elif gnew[ind] == 'NEG': 248 gnew[ind:ind+1] = [g for g in groups if g < 0] 249 else: 250 try: gnew[ind] = int(gnew[ind]) 251 except: 252 print("** invalid selection group '%s'" % gnew[ind]) 253 return 1 254 255 if self.verb > 2: print('-- red. by glist: groups %s' % gnew) 256 cols = self.ordered_cols_by_group_list(gnew) 257 if self.verb > 1: print('-- red. by glist: cols %s' % cols) 258 return self.reduce_by_vec_list(cols) 259 260 def show_header(self): 261 print('\n'.join(self.header)) 262 263 def show_group_labels(self): 264 show_groups = (len(self.groups) == self.nvec) 265 show_labs = (len(self.labels) == self.nvec) 266 if not show_groups and not show_labs: 267 print('** no label info to show') 268 return 269 270 for ind in range(self.nvec): 271 if self.verb: 272 if show_groups: gstr = ', group %-3s' % self.groups[ind] 273 else: gstr = '' 274 if show_labs: lstr = ', label %s' % self.labels[ind] 275 else: lstr = '' 276 print('index %3d%s%s' % (ind, gstr, lstr)) 277 elif show_labs: print('%s' % self.labels[ind]) 278 279 def show_df_info(self, protect=1): 280 """protect, just in case, since AP processing will depend on this""" 281 if self.verb > 4: print("++ show_df_info, protect=%d"%protect) 282 if protect: self.protected_show_df_info() 283 else: self._show_df_info() 284 285 def protected_show_df_info(self): 286 try: self._show_df_info() 287 except: print("** failed show_df_info") 288 289 def _show_df_info(self): 290 show_groups = (len(self.groups) == self.nvec) 291 show_labs = (len(self.labels) == self.nvec) 292 if not show_groups and not show_labs: 293 print('** no DF info to show') 294 return 295 if self.GLapplied != 1: 296 print("** show_df_info: requires original Xmat from 3dDeconvolve") 297 return 298 299 # known -ortvec label prefix list: ROI ROIPC bandpass morts mot 300 # seen by: 301 # grep -h ' -ortvec' results.*/proc.* | awk '{print $3}' \ 302 # | tr '._' ' ' | awk '{print $1}' | sort | uniq 303 304 305 ort_labs = ['ROIPC', 'ROI', 'bandpass', 'morts'] 306 mot_labs = ['mot_deriv', 'mot', 'roll', 'pitch', 'yaw', 'dS', 'dL', 'dP'] 307 nortlabs = len(ort_labs) 308 nmotlabs = len(mot_labs) 309 # put them together for extraction 310 all_ort_labs = ort_labs[:] 311 all_ort_labs.extend(mot_labs[:]) 312 313 # partition regressors by group (polort, other RONI, ROI) 314 pol_list = [ind for ind in range(self.nvec) if self.groups[ind] < 0] 315 zero_list = [ind for ind in range(self.nvec) if self.groups[ind] == 0] 316 pos_list = [ind for ind in range(self.nvec) if self.groups[ind] > 0] 317 318 bp_list = [ind for ind in zero_list \ 319 if UTIL.starts_with(self.labels[ind], "bandpass")] 320 321 # partition group 0 (RONI) 322 missed, ort_counts = self.count_sublist_labs(self.labels, all_ort_labs, 323 index_list=zero_list) 324 325 # get counts from any motion labels 326 # (now we have any ort_lab count, motcount, and missed) 327 motcount = UTIL.loc_sum(ort_counts[nortlabs:]) 328 329 # maybe Paul will want BP counts across runs 330 bp_run_counts = self.get_bp_counts_by_run(bp_list, verb=(self.verb>1)) 331 332 # make a print list of 3 element strings 333 plist = [] 334 335 # real info 336 ncensor = self.nrowfull - self.nt 337 totalDF = self.nrowfull 338 plist.append(self._df_entry("initial DF", totalDF, 1.0)) 339 plist.append([]) # blank line 340 # group > 0, censoring, -1 341 plist.append(self._df_entry("DF used for regs of interest", len(pos_list), 342 1.0*len(pos_list)/totalDF)) 343 plist.append(self._df_entry("DF used for censoring", ncensor, 344 1.0*ncensor/totalDF)) 345 plist.append(self._df_entry("DF used for polort", len(pol_list), 346 1.0*len(pol_list)/totalDF)) 347 # group 0: motion ort_lab entries, other (missed) 348 plist.append(self._df_entry("DF used for motion", motcount, 349 1.0*motcount/totalDF)) 350 351 # other counts are reported only if positive 352 for oind, olab in enumerate(ort_labs): 353 if ort_counts[oind] > 0: 354 plist.append(self._df_entry("DF used for %s"%olab, 355 ort_counts[oind], 1.0*ort_counts[oind]/totalDF)) 356 if missed > 0: 357 plist.append(self._df_entry("DF used for other RONI", missed, 358 1.0*missed/totalDF)) 359 360 plist.append(self._df_entry("total DF used", ncensor+self.nvec, 361 1.0*(ncensor+self.nvec)/totalDF)) 362 plist.append([]) # blank line 363 plist.append(self._df_entry("final DF", self.nt-self.nvec, 364 1.0*(self.nt-self.nvec)/totalDF)) 365 366 # get max string lengths for padding 367 m0 = max([len(p[0]) for p in plist if len(p) > 0]) 368 m1 = max([len(p[1]) for p in plist if len(p) > 0]) 369 m2 = max([len(p[2]) for p in plist if len(p) > 0]) 370 371 # and finally, print this out, assuming conversion to percent 372 print("") 373 for p in plist: 374 if len(p) == 3: 375 print("%-*s : %*s : %*s%%" % (m0, p[0], m1, p[1], m2, p[2])) 376 else: 377 print("") 378 print("") 379 380 def _df_entry(self, v1, v2, v3, v3_f2pct=1): 381 """return list of values converted to strings 382 if v3_f2pct: convert as fraction to percent 383 """ 384 if v3_f2pct: s3 = '%.1f'%(100.0*v3) 385 else: s3 = '%s' % v3 386 return ['%s'%v1, '%s'%v2, s3] 387 388 def get_bp_counts_by_run(self, bp_list, verb=0): 389 """return a list of bp_counts across runs""" 390 391 # note the index at the start of each run 392 run_starts = [0]*self.nruns 393 roff = 0 394 for rind, rlen in enumerate(self.run_len): 395 run_starts[rind] += roff 396 roff += rlen 397 if verb: print("run_starts : %s" % UTIL.int_list_string(run_starts)) 398 399 # for each bp reg, see which run it seems to belong to 400 bp_counts = [0]*self.nruns 401 for bind, bpind in enumerate(bp_list): 402 rind = self.get_bp_run(self.mat[bpind], run_starts) 403 if rind < 0 or rind >= self.nruns: 404 if verb: print("** bad bp rind = %d" % rind) 405 else: 406 bp_counts[rind] += 1 407 if verb: print("bp_counts : %s" % UTIL.int_list_string(bp_counts)) 408 409 return bp_counts 410 411 def count_sublist_labs(self, labels, prefix_list, index_list=None): 412 """given a list of labels and a corresponding index sub-list (None?), 413 return a list of counts of each prefix in prefix_list (checked in 414 order) 415 416 If index list is None, use all. 417 418 return num_not_found, and count list 419 """ 420 if not labels or len(labels) == 0: return 0, [] 421 if not prefix_list or len(prefix_list) == 0: return 0, [] 422 423 # first note which elements of 'labels' we will go after 424 if index_list == None: ind_list = list(range(len(labels))) 425 else: ind_list = index_list 426 427 notfound = 0 428 counts = [0]*len(prefix_list) 429 for lind in ind_list: 430 found = 0 431 for pind, prefix in enumerate(prefix_list): 432 if UTIL.starts_with(labels[lind], prefix): 433 counts[pind] += 1 434 found = 1 435 break 436 if not found: notfound += 1 437 438 return notfound, counts 439 440 def get_bp_run(self, bpreg, run_starts): 441 """given a bandpass regressor and the run starts, return the 0-based 442 run index that seems to apply 443 444 So find the min and max indices of non-zero value, and check that 445 they are in the same run. 446 """ 447 # for some reason, zeros are not exactly zero in the X-matrix 448 # (epsilon does not need to be too small, since all regressors 449 # should have a max close to 1) 450 epsilon = 0.001 451 nzmax = -1 452 453 nt = len(bpreg) 454 455 # find the first non-zero value 456 nzmin = -1 457 for tp in range(nt): 458 if abs(bpreg[tp]) > epsilon: 459 nzmin = tp 460 break 461 462 # and the last 463 nzmax = -1 464 for tp in range(nt-1, -1, -1): 465 if abs(bpreg[tp]) > epsilon: 466 nzmax = tp 467 break 468 469 # if we failed, bail 470 if nzmin < 0 or nzmax < 0: 471 print("** get_bp_run: BP reg is all zero?") 472 return -1 473 474 # see if they were in the same run 475 rind_min = self.get_run_index(run_starts, nzmin) 476 rind_max = self.get_run_index(run_starts, nzmax) 477 478 if rind_min != rind_max: 479 print("** get_bp_run: rind_min != rind_max (%d != %d)" \ 480 % (rind_min, rind_max)) 481 return -1 482 483 return rind_min 484 485 def get_run_index(self, run_starts, index): 486 """ignore first run start, and see if index is less than any""" 487 488 for rind in range(len(run_starts)-1): 489 if index < run_starts[rind+1]: 490 return rind 491 return len(run_starts) - 1 492 493 def reduce_by_label_prefix(self, keep_pre=[], drop_pre=[]): 494 495 rv, newlist = self.label_prefix_to_ints(keep_pre, drop_pre) 496 if rv: return 1 497 498 self.reduce_by_vec_list(newlist) 499 500 if self.verb>1: print('-- reduce_by_label_prefix, labels: %s'%self.labels) 501 502 return 0 503 504 def label_prefix_to_ints(self, keep_pre=[], drop_pre=[]): 505 """return a list of label indices, based on what is kept or dropped 506 keep labels in their original order 507 508 if keep_pre and one of keep_pre matches, keep 509 if drop_pre and one of drop_pre matches, remove 510 511 return 0 on success, along with the int list""" 512 513 if not self.ready: 514 print('** reduce_by_label_prefix: Afni1D is not ready') 515 return 1, [] 516 517 if len(self.labels) == 0: 518 print('** no dataset labels to reduce by...') 519 return 0, list(range(self.nvec)) 520 521 # make a list of label indices to keep, first, before removing 522 if len(keep_pre) > 0: 523 newlist = [] 524 for ind, lab in enumerate(self.labels): 525 for pre in keep_pre: 526 if UTIL.starts_with(lab, pre): 527 newlist.append(ind) 528 break 529 newlist = UTIL.get_unique_sublist(newlist) 530 if self.verb > 2: print('++ applied keep_pre to produce: %s' % newlist) 531 else: 532 newlist = list(range(len(self.labels))) 533 534 if len(drop_pre) > 0: 535 new2 = [] 536 for ind in newlist: 537 keep = 1 538 for pre in drop_pre: 539 if UTIL.starts_with(self.labels[ind], pre): 540 keep = 0 541 break 542 if keep: new2.append(ind) 543 newlist = new2 544 if self.verb > 2: print('++ applied drop_pre to produce %s' % newlist) 545 546 return 0, newlist 547 548 def run_info_is_consistent(self, whine=1): 549 """verify consistency of nruns/run_len/nt""" 550 if not self.ready: 551 if whine: print('** RIIC: not ready (to test consistency)') 552 return 0 553 if len(self.run_len) == 0: 554 if whine: print('** RIIC: ready be len(run_len) == 0!') 555 return 0 556 nt = UTIL.loc_sum(self.run_len) 557 if nt != self.nt: 558 if whine: 559 print('** RIIC: nt=%d != sum of run_len: %s'%(self.nt,self.run_len)) 560 return 0 561 return 1 562 563 # --- begin: mathematical manipulations --- 564 565 def transpose(self): 566 """transpose the matrix and swap nrows and ncols""" 567 if self.verb > 3: print('-- Afni1D transpose...') 568 if not self.ready: 569 print("** matrix '%s' is not ready for transposing" % self.name) 570 return 571 self.mat = [[row[i] for row in self.mat] for i in range(self.nt)] 572 newnt = self.nvec 573 self.nvec = self.nt 574 self.nt = newnt 575 self.run_len = [self.nt] # init to 1 run 576 577 def demean(self): 578 """demean each vector per run (mean of each run will be zero) 579 580 return 0 on success""" 581 582 if self.verb > 3: print('-- Afni1D demean...') 583 584 if not self.ready: 585 print('** demean: Afni1D is not ready') 586 return 1 587 588 # verify nruns and run_len 589 if not self.run_info_is_consistent(whine=1): 590 if self.verb > 1: print('** runs inconsistent for demean') 591 return 1 592 593 if self.verb > 1: 594 print("-- demean: over %d runs of lengths %s" \ 595 % (self.nruns, self.run_len)) 596 597 # foreach run of each vector, demean 598 for ind in range(self.nvec): 599 offset = 0 600 for t in self.run_len: 601 UTIL.demean(self.mat[ind], offset, offset+t-1) 602 offset += t 603 604 return 0 605 606 def derivative(self, direct=0): 607 """change each value to its derivative (cur val - previous) 608 new time[0] will always be 0 609 610 process the data per run, so derivatives do not cross boundaries 611 612 direct: 0 means backward difference, 1 means forward 613 614 return 0 on success""" 615 616 if self.verb > 3: print('-- Afni1D derivative...') 617 618 if not self.ready: 619 print('** derivative: Afni1D is not ready') 620 return 1 621 622 # verify nruns and run_len 623 if not self.run_info_is_consistent(whine=1): 624 if self.verb > 1: print('** runs inconsistent for derivative') 625 return 1 626 627 if self.verb > 1: 628 print("-- derivative: over %d runs of lengths %s" \ 629 % (self.nruns, self.run_len)) 630 631 # apply derivative to each run of each vector 632 for ind in range(self.nvec): 633 UTIL.derivative(self.mat[ind], in_place=1, direct=direct) 634 # if forward diff, clear final indexes, else clear initial 635 if direct: 636 offset = -1 637 for t in self.run_len: 638 offset += t 639 self.mat[ind][offset] = 0 640 else: 641 offset = 0 642 for t in self.run_len: 643 self.mat[ind][offset] = 0 644 offset += t 645 646 return 0 647 648 def abs(self): 649 """take the absolute value of each entry 650 return 0 on success, 1 on error""" 651 652 if self.verb > 3: print('-- Afni1D abs...') 653 654 if not self.ready: 655 print('** abs: Afni1D is not ready') 656 return 1 657 658 for ind in range(self.nvec): 659 newvec = [] 660 for rind in range(self.nt): 661 self.mat[ind][rind] = abs(self.mat[ind][rind]) 662 663 return 0 664 665 def bool_negate(self): 666 """set clear values and clear set values 667 return 0 on success, 1 on error""" 668 669 if self.verb > 3: print('-- bool_negate: old zeros = %d' \ 670 % self.mat[0].count(0)) 671 672 if not self.ready: 673 print('** bool_negate: Afni1D is not ready') 674 return 1 675 676 for v in range(self.nvec): 677 for t in range(self.nt): 678 if self.mat[v][t]: self.mat[v][t] = 0 679 else: self.mat[v][t] = 1 680 681 if self.verb > 3: print('-- negate: new zeros = %d' % self.mat[0].count(0)) 682 683 return 0 684 685 def volreg_2_allineate(self): 686 """convert a 6-parameter time series from 3dvolreg to a 12-parameter 687 time series for 3dAllineate 688 689 params (3dvolreg -1Dfile): roll, pitch, yaw, dS, dL, dP 690 (3dAllineate 1Dapply): 691 692 3dvolreg params: v0 v1 v2 v3 v4 v5 693 3dAllieate params: -v4 -v5 -v3 v0 v1 v2 694 695 Permute and negate vectors, and append [0] vectors. 696 697 return 0 on success 698 """ 699 if self.verb > 3: print('-- volreg_2_allineate...') 700 if not self.ready: 701 print("** matrix '%s' is not ready for volreg2allin" % self.name) 702 return 1 703 if self.nvec != 6: 704 print("** matrix '%s' does not have 6 columns" % self.name) 705 return 1 706 707 mm = self.mat # for convenience 708 709 # negate second triplet of vectors 710 for col in range(3,6): 711 mm[col] = [-val for val in mm[col]] 712 713 # permute and append 6 [0] vectors 714 zvec = [0] * self.nt 715 self.mat = [mm[4], mm[5], mm[3], mm[0], mm[1], mm[2], 716 zvec, zvec, zvec, zvec, zvec, zvec ] 717 self.nvec = 12 718 719 return 0 720 721 def collapse_cols(self, method, weight=None): 722 """collapsed the matrix to a single array of length nt (nvec will = 1) 723 724 collapsing will apply 'method' across the nvec per time point 725 726 method = 'min' : min across nvec 727 method = 'minabs' : min abs across nvec 728 method = 'max' : max across nvec 729 method = 'maxabs' : max abs across nvec 730 method = 'euclidean_norm' : sqrt(sum squares) 731 or = 'enorm' 732 method = 'weighted_enorm' : sqrt(sum weighted squares) 733 734 Note: the result will still be a trivial 2-D array, where element 0 735 is the collapsed time series. This allows other functionality 736 to still work. 737 738 return 0 on success, 1 on error""" 739 740 if self.verb > 3: print('-- collapse cols, method = %s' % method) 741 742 if not self.ready: 743 print('** collapse: Afni1D is not ready') 744 return 1 745 746 # either create a new mat and apply later, or modify mat and recur 747 748 if method == 'min': 749 mat = [min([self.mat[v][t] for v in range(self.nvec)]) \ 750 for t in range(self.nt)] 751 elif method == 'max': 752 mat = [max([self.mat[v][t] for v in range(self.nvec)]) \ 753 for t in range(self.nt)] 754 elif method == 'minabs': # take abs and recur 755 if self.abs(): return 1 756 return self.collapse_cols('min') 757 elif method == 'maxabs': # take abs and recur 758 if self.abs(): return 1 759 return self.collapse_cols('max') 760 elif method == 'enorm' or method == 'euclidean_norm': 761 mat = [UTIL.euclidean_norm([self.mat[v][t] for v in range(self.nvec)])\ 762 for t in range(self.nt)] 763 elif method == 'weighted_enorm': 764 if weight == None: 765 print('** missing weight vector for weighted_enorm') 766 return 1 767 if len(weight) != self.nvec: 768 print('** weighted_enorm weight vector length (%d) != nvec (%d)' \ 769 % (len(weight), self.nvec)) 770 return 1 771 mat = [UTIL.weighted_enorm( 772 [self.mat[v][t] for v in range(self.nvec)], weight) \ 773 for t in range(self.nt)] 774 else: 775 print("** collapse_cols: unknown method:", method) 776 return 1 777 778 # note that there is only a single column left 779 self.nvec = 1 780 self.clear_cormat() 781 782 del(self.mat) 783 self.mat = [mat] 784 return 0 785 786 def mat_times_vec(self, vec): 787 """modify matrix as self.mat * vec""" 788 lv = len(vec) 789 if lv == 0 or self.nt == 0 or self.nvec == 0: 790 if self.verb > 1: print('** have mat_times_vec with empty inputs') 791 return 792 793 for col, mvec in enumerate(self.mat): 794 val = vec[col] 795 for ind in range(self.nt): 796 mvec[ind] *= val 797 798 def unitize(self): 799 """make every vector unit length, use existing memory""" 800 801 for ind, vec in enumerate(self.mat): 802 vlen = UTIL.L2_norm(vec) 803 if vlen == 0.0: 804 for ind in range(self.nt): vec[ind] = 0 805 else: 806 for ind in range(self.nt): vec[ind] /= vlen 807 808 def project_out_vec(self, vec=None, unit=0): 809 """project a vector out of the matrix vectors, if None, use mean 810 811 for each vector y, subtract <vT,y>/<vT.v> * v 812 813 if unit: subtract mean of unit vectors 814 815 return 0 on success, else 1 816 """ 817 if vec == None: 818 if unit: vec = self.get_mean_unit_vec() 819 else: vec = self.get_mean_vec() 820 if len(vec) != self.nt: 821 print('** project_out_vec: bad length %d != %d' % (len(vec), self.nt)) 822 return 1 823 824 vlen = UTIL.L2_norm(vec) 825 if vlen == 0: # then nothing to do 826 if self.verb > 0: print('-- project_out_vec: empty vector to remove') 827 return 0 828 vunit = UTIL.lin_vec_sum(1.0/vlen, vec, 0, None) 829 830 self.mat = [UTIL.proj_out_vec(mv, vunit, unit_v2=1) for mv in self.mat] 831 832 # --- end: mathematical manipulations --- 833 834 # --- these functions return some aspect of the matrix --- 835 836 def get_mean_vec(self): 837 """return a single vector as the mean of all matrix vectors""" 838 v = [0] * self.nt 839 if self.nt == 0: return [] 840 for ind in range(self.nt): 841 v[ind] = UTIL.loc_sum([self.mat[i][ind] for i in range(self.nvec)]) 842 return UTIL.lin_vec_sum(1.0/self.nvec, v, 0, None) 843 844 def get_mean_unit_vec(self, demean=1): 845 """return a single vector as the mean of all matrix vectors 846 AFTER the vectors have been (demeaned? and) scaled to unit lengths 847 """ 848 adcopy = self.copy() 849 if demean: adcopy.demean() 850 adcopy.unitize() 851 gu = adcopy.get_mean_vec() 852 del(adcopy) 853 return gu 854 855 def show_gcor_all(self): 856 for ind in range(1,6): 857 print('----------------- GCOR test %d --------------------' % ind) 858 exec('val = self.gcor%d()' % ind) 859 exec('print("GCOR rv = %.10f" % val)') 860 861 def show_gcor_doc_all(self): 862 for ind in range(1,6): 863 print('----------------- GCOR doc %d --------------------' % ind) 864 exec('val = self.gcor%d.__doc__' % ind) 865 exec('print("%s" % val)') 866 867 # basically, a link to the one we really want to call 868 def gcor(self): return self.gcor2() 869 870 def show_gcor(self, verb=-1): 871 gc = self.gcor() 872 nv = self.nvec 873 nt = self.nt 874 if verb < 0: verb = self.verb 875 if verb: print("GCOR for %d time series of length %d = %g" % (nv, nt, gc)) 876 else: print("%g" % gc) 877 878 def gcor1(self): 879 """GOLD STANDARD 880 881 compute GCOR via: 882 - fill cormat (correlation matrix) 883 - return average 884 """ 885 self.set_cormat() 886 if not self.cormat_ready: return 0 887 888 ss = 0.0 889 for r in range(self.nvec): 890 for c in range(self.nvec): 891 ss += self.cormat[r][c] 892 893 return float(ss)/(self.nvec*self.nvec) 894 895 def gcor2(self): 896 """PLATNUM STANDARD (method applied in afni_proc.py) 897 898 compute GCOR via: 899 - demean 900 - unitize 901 - get average unit time series gu 902 - return sum_of_squares(gu) = length(gu)^2 903 """ 904 adcopy = self.copy() 905 adcopy.demean() 906 adcopy.unitize() 907 gu = adcopy.get_mean_vec() 908 909 en = UTIL.dotprod(gu,gu) 910 911 return en 912 913 def gcor3(self): 914 """compute GCOR via: 915 - demean 916 - unitize 917 - get average unit time series gu 918 - compute average correlation of gu and each unit vector 919 - return |gu| times ave corr with gu 920 """ 921 922 adcopy = self.copy() 923 adcopy.demean() 924 adcopy.unitize() 925 926 gu = copy.deepcopy(adcopy.mat[0]) 927 for ind in range(1, adcopy.nvec): 928 gu = UTIL.lin_vec_sum(1, gu, 1, adcopy.mat[ind]) 929 gu = [val/float(adcopy.nvec) for val in gu] 930 931 ss = 0.0 932 for r in range(adcopy.nvec): 933 ss += UTIL.correlation_p(gu, adcopy.mat[r]) 934 935 ss /= adcopy.nvec 936 937 return UTIL.euclidean_norm(gu)*ss 938 939 def gcor4(self): 940 """compute GMEAN via: 941 - get average time series gmean 942 - get average correlation of gmean and each unit vector 943 - return square (akin to above, but with gmean, not gu) 944 """ 945 946 ss = self.get_ave_correlation_w_vec(self.get_mean_vec()) 947 948 print("ave corr = %.10f, square = %.10f" % (ss, ss*ss)) 949 950 return ss*ss 951 952 def gcor5(self): 953 """compute GCOR via: 954 - get average time series mean, gmean 955 - get average unit time series mean, gu 956 - return (gu.gmean/|gmean|)^2 957 """ 958 959 adcopy = self.copy() 960 adcopy.demean() 961 962 m0 = adcopy.get_mean_vec() 963 964 adcopy.unitize() 965 m1 = adcopy.get_mean_vec() 966 967 dd = UTIL.dotprod(m0, m1) 968 l0 = UTIL.euclidean_norm(m0) 969 l1 = UTIL.euclidean_norm(m1) 970 c = UTIL.correlation_p(m0, m1) 971 972 print("len(gmean) = %.10f, square = %.10f" % (l0, l0*l0)) 973 print("len(gunit) = %.10f, square = %.10f" % (l1, l1*l1)) 974 print("corr(gm, gu) = %.10f" % (c)) 975 print("corr(gm, gu)*len(gu) = %.10f" % (c*l1)) 976 print("squared = %.10f" % (c*c*l1*l1)) 977 978 return dd*dd/(l0*l0) 979 980 def get_ave_correlation_w_vec(self, vec): 981 """return the average correlation of each vector with vec""" 982 gu = self.get_mean_vec() 983 ss = 0 984 for r in range(self.nvec): 985 ss += UTIL.correlation_p(vec, self.mat[r]) 986 ss /= self.nvec 987 988 return ss 989 990 def get_gcor_wout_gmean(self, unitize=0): 991 """project out global mean, then compute gcor 992 unitize: unitize before projection 993 """ 994 adcopy = self.copy() 995 if unitize: adcopy.unitize() 996 adcopy.project_out_vec() 997 gc = adcopy.gcor() 998 del(adcopy) 999 return gc 1000 1001 # --- end: functions returning some aspect of the matrix --- 1002 1003 def copy(self): 1004 """return a complete (deep)copy of the current object""" 1005 return copy.deepcopy(self) 1006 1007 def extreme_mask(self, emin, emax, inclusive=1): 1008 """convert to a time series mask of 0/1 values where the result is 1009 1, if value is extreme (>emax or <emin, with = if inclusive) 1010 0, otherwise (moderate) 1011 1012 For example, these would be the TRs to omit in a censor.1D file. 1013 1014 return 0 on success""" 1015 1016 if self.verb > 3: 1017 if inclusive: print('-- extreme_mask: excluding (%g,%g)'%(emin,emax)) 1018 else: print('-- extreme_mask: excluding [%g,%g]'%(emin,emax)) 1019 1020 if not self.ready: 1021 print('** extreme_mask: Afni1D is not ready') 1022 return 1 1023 1024 if emin > emax: 1025 print('** extreme_mask: emin > emax (', emin, emax, ')') 1026 return 1 1027 1028 self.mat = [UTIL.vec_extremes(vec,emin,emax,inclusive)[1] \ 1029 for vec in self.mat] 1030 1031 if self.verb > 1: 1032 count = UTIL.loc_sum([val for val in self.mat[0] if val == 1]) 1033 print('++ extreme_mask: removing %d of %d vals' % (count,self.nt)) 1034 1035 return 0 1036 1037 def moderate_mask(self, emin, emax, inclusive=1): 1038 """convert to a time series mask of 0/1 values where the result is 1039 1, if value is moderate (within (emin,emask), with = if inclusive) 1040 0, otherwise (extreme) 1041 1042 For example, these would be the TRs to keep in a censor.1D file. 1043 1044 return 0 on success""" 1045 1046 if self.verb > 3: 1047 if inclusive: print('-- moderate_mask: keeping (%g,%g)'%(emin,emax)) 1048 else: print('-- moderate_mask: keeping [%g,%g]'%(emin,emax)) 1049 1050 if not self.ready: 1051 print('** moderate_mask: Afni1D is not ready') 1052 return 1 1053 1054 if emin > emax: 1055 print('** moderate_mask: emin > emax (', emin, emax, ')') 1056 return 1 1057 1058 self.mat = [UTIL.vec_moderates(vec,emin,emax,inclusive)[1] \ 1059 for vec in self.mat] 1060 1061 if self.verb > 1: 1062 count = UTIL.loc_sum([val for val in self.mat[0] if val == 1]) 1063 print('++ moderate_mask: keeping %d of %d vals' % (count,self.nt)) 1064 1065 return 0 1066 1067 def set_first_TRs(self, nfirst, newval=1): 1068 """set the first nfirst TRs to newval""" 1069 1070 if self.verb > 3: 1071 print('-- setting first %d TRs, newval = %d ...' % (nfirst, newval)) 1072 1073 if not self.ready: 1074 print('** set_first_TRs: Afni1D is not ready') 1075 return 1 1076 1077 # apply 'newval' to first 'nfirst' in each run 1078 for ind in range(self.nvec): 1079 offset = 0 1080 for run, rlen in enumerate(self.run_len): 1081 if nfirst > rlen: 1082 print('** mask_first_TRs, nfirst %d > run len %d (of run %d)' \ 1083 % (nfirst, rlen, run+1)) 1084 return 1 1085 # apply censor val for first nfirst TRs 1086 for tr in range(nfirst): self.mat[ind][offset+tr] = newval 1087 offset += rlen 1088 1089 return 0 1090 1091 def mask_prior_TRs(self): 1092 """if one TR is set, also set the prior one""" 1093 1094 if self.verb > 3: print('-- masking prior TRs...') 1095 1096 if not self.ready: 1097 print('** mask_prior_TRs: Afni1D is not ready') 1098 return 1 1099 1100 for v in range(self.nvec): 1101 for t in range(self.nt-1): 1102 if self.mat[v][t+1] and not self.mat[v][t]: self.mat[v][t] = 1 1103 1104 return 0 1105 1106 def clear_next_TRs(self): 1107 """if one TR is clear, also clear the next one""" 1108 1109 if self.verb > 3: print('-- clearing next TRs...') 1110 1111 if not self.ready: 1112 print('** clear_next_TRs: Afni1D is not ready') 1113 return 1 1114 1115 for v in range(self.nvec): 1116 for t in range(self.nt-2, -1, -1): 1117 if not self.mat[v][t] and self.mat[v][t+1]: self.mat[v][t+1] = 0 1118 1119 return 0 1120 1121 def clear_prior_TRs(self): 1122 """if one TR is clear, also clear the prior one""" 1123 1124 if self.verb > 3: print('-- clearing prior TRs...') 1125 1126 if not self.ready: 1127 print('** clear_prior_TRs: Afni1D is not ready') 1128 return 1 1129 1130 for v in range(self.nvec): 1131 for t in range(self.nt-1): 1132 if not self.mat[v][t+1] and self.mat[v][t]: self.mat[v][t] = 0 1133 1134 return 0 1135 1136 def pad_into_many_runs(self, rindex1, nruns, rlengths=[]): 1137 """pad over time so that this is run #rindex1 out of nruns runs 1138 (rindex1 should be 1-based) 1139 1140 if rlengths is not set, assume all runs are of a set length 1141 if rlengths IS set, nruns is ignored 1142 1143 Only one of nruns/rlengths is is applied. 1144 1145 return 0 on success""" 1146 1147 if self.verb > 3: print('-- pad_into_many_runs: rindex1=%d, nruns=%d' \ 1148 % (rindex1,nruns)) 1149 1150 if not self.ready: 1151 print('** pad into runs: Afni1D is not ready') 1152 return 1 1153 1154 # ---- set NR (nruns) and rlens (list of lengths) 1155 1156 # decide whether to use rlengths or nruns 1157 if len(rlengths) > 0: 1158 NR = len(rlengths) 1159 rlens = rlengths 1160 else: # assume nruns runs of length self.nt 1161 if nruns < 1: 1162 print('** pad into runs: bad nruns (%d)' % nruns) 1163 return 1 1164 NR = nruns 1165 rlens = [self.nt for r in range(NR)] 1166 1167 # ---- check for consistency 1168 1169 # verify rindex1 (using 1-based run index) 1170 if rindex1 < 1 or rindex1 > NR: 1171 print('** pad into runs: run index (%d) out of range [1,%d]' \ 1172 % (rindex1, NR)) 1173 return 1 1174 1175 # apply rind as 0-based run index 1176 rind = rindex1-1 1177 1178 # if rlengths, verify match for run #rindex1 1179 if self.nt != rlens[rind]: 1180 print("** cannot pad into many runs, nt (%d) != rlens[%d] (%d)" \ 1181 % (self.nt, rind, rlens[rind])) 1182 return 1 1183 1184 # ---- and do the work 1185 1186 # first just stick in a 0 at time t=0 1187 for row in range(self.nvec): 1188 self.mat[row] = [0] + self.mat[row] 1189 self.nt += 1 1190 1191 # now create a time list that sticks t0 in as padding 1192 # (fill runs before rind, then rind, then after) 1193 tlist = [] 1194 for run in range(rind): 1195 tlist.extend([0 for t in range(rlens[run])]) 1196 tlist.extend([t+1 for t in range(rlens[rind])]) # skipping pad at t=0 1197 for run in range(NR-1-rind): 1198 tlist.extend([0 for t in range(rlens[run+rind+1])]) 1199 1200 # and apply tlist 1201 if self.reduce_by_tlist(tlist): return 1 1202 1203 # update run info 1204 self.nruns = NR 1205 self.run_len = rlens 1206 self.nt = UTIL.loc_sum(rlens) 1207 1208 return 0 1209 1210 def randomize_trs(self, seed=0): 1211 """reorder the matrix columns randomly""" 1212 if self.verb > 1: print('-- randomizing %d trs (seed=%d)' \ 1213 % (self.nt, seed)) 1214 if seed: 1215 import random 1216 random.seed(seed) 1217 for vec in self.mat: UTIL.shuffle(vec) 1218 1219 def sort(self, reverse=0): 1220 """sort data over time axis (possibly reverse order)""" 1221 1222 if self.verb > 3: print('-- Afni1D sorting...') 1223 1224 for ind in range(self.nvec): 1225 self.mat[ind].sort(reverse=reverse) 1226 1227 return 0 1228 1229 def reverse(self): 1230 """reverse data over time axis""" 1231 1232 if self.verb > 3: print('-- Afni1D reverse...') 1233 1234 ilist = UTIL.decode_1D_ints('$..0(-1)',verb=self.verb,imax=self.nt-1) 1235 if self.reduce_by_tlist(ilist): return 1 1236 1237 return 0 1238 1239 def get_allzero_cols(self): 1240 """return a list of column indices for which the column is all zero 1241 (these might represent regressors which were censored out) 1242 1243 return status (0=success) and the column index list 1244 """ 1245 1246 return 0, [ind for ind in range(self.nvec) \ 1247 if UTIL.vals_are_constant(self.mat[ind], 0)] 1248 1249 return 0, rlist 1250 1251 def get_censored_trs(self, run_index=-1): 1252 """return a list of TRs that were censored 1253 (basically, return an inverted goodlist) 1254 1255 if run_index >= 0: restrict to that 0-based run 1256 (apply runstart info after invert_int_list) 1257 1258 return status (0=success) and the TR index list 1259 """ 1260 1261 rv, ulist = self.get_uncensored_trs() # no run_index 1262 if rv: return 1, [] 1263 1264 if self.nrowfull > 0: nt = self.nrowfull 1265 else: nt = self.nt 1266 1267 tr_list = UTIL.invert_int_list(ulist, top=nt-1) 1268 1269 return self.restrict_tr_list_to_run(tr_list, run_index) 1270 1271 def get_uncensored_trs(self, run_index=-1): 1272 """return a list of TRs that were used, i.e. we not censored 1273 (basically, return goodlist) 1274 1275 there are 2 valid types of inputs: 1276 1277 1. X-matrix format 1278 len(goodlist) > 0 and nrowfull >= len(goodlist) 1279 1280 if run_index >= 0, limit goodlist indices base on run_len[] 1281 (run index is 0-based) 1282 1283 2. binary format 1284 nvec == 1 and UTIL.vals_are_0_1(mat[0]) 1285 1286 return status (0=success) and the TR index list 1287 """ 1288 1289 if not self.ready: 1290 print("** Afni1D not ready for get_uncensored_trs") 1291 return 1, [] 1292 1293 # either have goodlist from xmat or take indices of non-zero values 1294 if len(self.goodlist) > 0: 1295 tr_list = self.goodlist 1296 else: 1297 tr_list = [i for i,v in enumerate(self.mat[0]) if v] 1298 1299 # and restrict to run_index 1300 return self.restrict_tr_list_to_run(tr_list, run_index) 1301 1302 def uncensor_from_vector(self, cvec): 1303 """zero-pad list to match length of cvec 1304 current list length must match uncensored cved list 1305 """ 1306 if self.nruns > 1: 1307 print('** cannot uncensor multiple runs from vector') 1308 return 1 1309 1310 clen = len(cvec) 1311 ncen = cvec.count(0) 1312 if self.nt != clen - ncen: 1313 print('** uncensor from vec: nt = %d, but nocen len = %d' \ 1314 % (self.nt, clen-ncen)) 1315 return 1 1316 1317 # explicitly set newvec, including with zeros 1318 newmat = [] 1319 for mvec in self.mat: 1320 newvec = cvec[:] 1321 goodind = 0 1322 for ind, val in enumerate(cvec): 1323 if val: 1324 newvec[ind] = mvec[goodind] 1325 goodind += 1 1326 else: 1327 newvec[ind] = 0 1328 newmat.append(newvec) 1329 1330 if self.verb > 1: 1331 print('++ uncensored %d values from vector' % ncen) 1332 if self.verb > 2: 1333 cenlist = ['%s'%i for i in range(clen) if not cvec[i]] 1334 print(' zero-pad list: %s' % (', '.join(cenlist))) 1335 del(cenlist) 1336 1337 del(self.mat) 1338 self.mat = newmat 1339 self.nt = clen 1340 self.run_len = [self.nt] 1341 1342 return 0 1343 1344 def restrict_tr_list_to_run(self, tr_list, run_index=-1): 1345 """rely on runstart to restrict tr_list 1346 return status and restricted list 1347 """ 1348 1349 # maybe there is nothing to do 1350 if run_index < 0: return 0, tr_list 1351 1352 if run_index >= self.nruns: 1353 print('** cannot restrict TR list for (1-based) run %d, have %d' \ 1354 % (run_index+1, self.nruns)) 1355 return 1, [] 1356 1357 if self.nrowfull > 0: nt = self.nrowfull 1358 else: nt = self.nt 1359 1360 # restrict return values to [bot, top) 1361 bot = self.runstart[run_index] 1362 if run_index == self.nruns - 1: top = nt 1363 else: top = self.runstart[run_index+1] 1364 1365 if self.verb > 2: 1366 print('-- restricting %d TRs between %d and %d' \ 1367 % (len(tr_list),bot,top-1)) 1368 1369 return 0, [v for v in tr_list if v >= bot and v < top] 1370 1371 def show_censor_count(self, invert=0, column=0): 1372 """display the total number of TRs censored (clear) in the given column 1373 1374 If multi-column data, can choose one. 1375 1376 invert : invert column count 1377 column : which column to count values from 1378 1379 return status""" 1380 1381 if self.verb > 3: 1382 print('++ showing censor count (inv=%d, col=%d)' % (invert, column)) 1383 1384 if not self.ready: 1385 print("** Afni1D not ready for write_timing to '%s'" % fname) 1386 return 1 1387 1388 ccount = self.mat[0].count(0) # start by counting 0 1389 lstr = '(simple form)' 1390 1391 # if it is a valid x-mat, override (should be a separate function) 1392 if self.havelabs: 1393 lstr = '(from xmat header)' 1394 rv, tplist = self.get_censored_trs() 1395 if not rv: 1396 ccount = len(tplist) 1397 1398 if invert: ccount = self.nt - ccount 1399 1400 if self.verb: print('total number of censored TRs %s = %d'%(lstr, ccount)) 1401 else: print(ccount) 1402 1403 return 0 1404 1405 # ---------------------------------------------------------------------- 1406 # cluster table functions 1407 1408 def is_csim_type(self): 1409 """test whether this element seems to be a clustsim table""" 1410 if len(self.header) < 1: return 0 1411 if len(self.mat) <= 1: return 0 1412 1413 if self.csim_fill_obj(): return 0 # if csim, status -> 1 1414 if self.csimstat < 1: return 0 1415 1416 return 1 1417 1418 def csim_show_min_clust_size(self, pthr=0.001, alpha=0.05, verb=1): 1419 """given pthr/alpha, show min cluster size and cluster details, 1420 depending on verb 1421 1422 verb 0: just print size 1423 1: also so attributes 1424 2: all debug attributes 1425 """ 1426 1427 csize = self.csim_clust_size(pthr, alpha) 1428 if verb == 0: 1429 print("%s" % csize) 1430 return 0 1431 1432 # else, show more verbose output 1433 self.csim_show_attrs(verb=verb) 1434 print(" pthr : %s" % pthr) 1435 print(" alpha : %s" % alpha) 1436 print(" min clust nvox : %s" % csize) 1437 1438 return 0 1439 1440 def csim_clust_size(self, pthr=0.001, alpha=0.05): 1441 """given pthr/alpha, return min clust size 1442 return 0 on error""" 1443 1444 if self.csim_fill_obj(): return 0 1445 if not self.csim_has_all_attrs(): return 0 1446 1447 # get to work 1448 try: 1449 pind = self.mat[0].index(pthr) 1450 except: 1451 print('** uncorrected pthr = %s not in ClustSim table' % pthr) 1452 return 0 1453 1454 try: 1455 aind = self.csimobj.avals.index(alpha) 1456 except: 1457 print('** corrected alpha = %s not in ClustSim table' % alpha) 1458 return 0 1459 1460 try: 1461 csize = int(math.ceil(self.mat[aind+1][pind])) 1462 except: 1463 print("** csim_clust_size: unknown index error") 1464 return 0 1465 1466 return csize 1467 1468 def csim_show_attrs(self, verb=1): 1469 if not self.csim_has_vo(): 1470 print("** ClustSim: no attribute object") 1471 return 1 1472 1473 if not self.csim_has_all_attrs(verb=verb): 1474 print("** ClustSim attributes are missing\n") 1475 1476 print("") 1477 1478 if verb > 1: 1479 self.csimobj.show() 1480 1481 # make a blur parameter estimate string 1482 param_str = '' 1483 try: 1484 params = self.csimobj.val('bvals') 1485 param_str = ', '.join(['%s' % p for p in params]) 1486 except: 1487 pass 1488 1489 print("ClustSim attributes:") 1490 print(" neighbors : NN-%s" % self.csimobj.val('NN')) 1491 print(" sidedness : %s" % self.csimobj.val('sided')) 1492 print(" blur est type : %s" % self.csimobj.val('btype')) 1493 print(" blur params : %s" % param_str) 1494 print(" grid voxels : %s" % self.csimobj.val('grid_nvox')) 1495 print(" grid vox size : %s" % self.csimobj.val('grid_vsize')) 1496 print(" mask : %s" % self.csimobj.val('mask')) 1497 print(" mask N voxels : %s" % self.csimobj.val('mask_nvox')) 1498 print("") 1499 1500 return 0 1501 1502 def csim_has_vo(self): 1503 if self.csimstat != 1: return 0 1504 if self.VO == None: return 0 1505 if not isinstance(self.csimobj, self.VO.VarsObject): return 0 1506 return 1 1507 1508 def csim_has_all_attrs(self, verb=1): 1509 if not self.csim_has_vo(): return 0 1510 1511 alist = ['command', 'btype', 'bvals', 'sided', 'grid_nvox', 'grid_vsize', 1512 'mask', 'mask_nvox', 'NN', 'avals'] 1513 cobj = self.csimobj 1514 attrs = cobj.attributes() 1515 hasall = 1 1516 for aname in alist: 1517 if not aname in attrs: 1518 hasall = 0 1519 if verb and not cobj.whined: 1520 print('** csim obj, missing attribute %s' % aname) 1521 1522 # only whine once 1523 if verb and not cobj.whined and not hasall: 1524 cobj.whined = 1 1525 1526 return hasall 1527 1528 # if we like this, formalize via class ClusterTable, say 1529 def csim_fill_obj(self, refill=0): 1530 # if not refill and we have set status, there is nothing to do 1531 if not refill and self.csimstat != -1: 1532 return 0 1533 1534 # fill as if we have not been here before 1535 self.csimstat = 0 1536 1537 if len(self.header) < 1: return 0 1538 if len(self.mat) <= 1: return 0 1539 if not UTIL.starts_with(self.header[0], '# 3dClustSim '): return 0 1540 1541 try: 1542 # try this out, it is separate, so make_random_timing need not import 1543 from afnipy import lib_vars_object as VO 1544 self.VO = VO 1545 cobj = self.VO.VarsObject() 1546 cobj.whined = 0 1547 for cline in self.header: 1548 csplit = cline.split() 1549 if UTIL.starts_with(cline, '# 3dClustSim '): 1550 # command, btype, bvals[3] 1551 cobj.command = cline[2:] 1552 for ind in range(len(csplit)): 1553 # get blur option 1554 if csplit[ind] == '-acf': 1555 cobj.btype = 'ACF' 1556 cobj.bvals = [float(v) for v in csplit[ind+1:ind+4]] 1557 elif csplit[ind] == '-fwhm': 1558 cobj.btype = 'FWHM' 1559 cobj.bvals = [float(v) for v in csplit[ind+1:ind+4]] 1560 elif csplit[ind] == '-fwhmxyz': 1561 cobj.btype = 'acf' 1562 bval = float(csplit[ind+1]) 1563 cobj.bvals = [bval, bval, bval] 1564 elif csplit[ind] == '-insdat': 1565 # if no mask, use this one 1566 if not cobj.valid('mask'): 1567 cobj.mask = csplit[ind+1] 1568 cobj.insdat = csplit[ind+1] 1569 cobj.btype = 'NONE' 1570 cobj.bvals = [] 1571 # 2 ways to get mask: -mask or -insdat 1572 elif csplit[ind] == '-mask': 1573 cobj.mask = csplit[ind+1] 1574 else: 1575 continue 1576 1577 elif cline.find('thresholding') > 1: 1578 cobj.sided = csplit[1] 1579 1580 elif UTIL.starts_with(cline, '# Grid: '): 1581 # grid_nvox, grid_vsize, mask_nvox 1582 gline = cline[8:] 1583 vind = gline.find(' (') 1584 gvals = gline[0:vind].split() 1585 cobj.grid_nvox = gvals[0] 1586 cobj.grid_vsize = ' '.join(gvals[1:3]) 1587 msplit = gline[vind+2:].split() 1588 cobj.mask_nvox = int(msplit[0]) 1589 1590 elif UTIL.starts_with(cline, '# -NN '): 1591 cobj.NN = int(csplit[2]) 1592 1593 elif UTIL.starts_with(cline, '# pthr '): 1594 cobj.avals = [float(csplit[i]) for i in range(3,len(csplit))] 1595 self.csimobj = cobj 1596 except: 1597 print('** failed to convert 3dClustSim header to object') 1598 return 1 1599 1600 self.csimstat = 1 1601 return 0 1602 1603 # ---------------------------------------------------------------------- 1604 1605 def write(self, fname, sep=" ", overwrite=0, with_header=0): 1606 """write the data to a new .1D file 1607 1608 fname : filename is required 1609 sep : separator between elements 1610 overwrite : whether to allow overwrite 1611 with_header : include comment-style header info 1612 1613 return status""" 1614 1615 if self.verb > 2: print('-- Afni1D write to %s, o=%s, h=%s' \ 1616 %(fname, overwrite, with_header)) 1617 1618 if not self.ready: 1619 print("** Afni1D not ready for write to '%s'" % fname) 1620 return 1 1621 if not fname: 1622 print("** missing filename for write") 1623 return 1 1624 1625 if fname == '-' or fname == 'stdout': fp = sys.stdout 1626 else: 1627 # normal file name: check existence and open 1628 if os.path.exists(fname) and not overwrite: 1629 print("** output file '%s' exists and 'overwrite' not set..."%fname) 1630 return 1 1631 1632 fp = open(fname, 'w') 1633 if not fp: 1634 print("** failed to open '%s' for writing" % fname) 1635 return 1 1636 1637 if with_header: 1638 hstr = self.make_xmat_header_string() 1639 if hstr != '': fp.write(hstr+'\n#\n') 1640 1641 for row in range(self.nt): 1642 fp.write(sep.join(['%g' % self.mat[i][row] for i in range(self.nvec)])) 1643 fp.write('\n') 1644 1645 fp.close() 1646 1647 return 0 1648 1649 def make_xmat_header_string(self): 1650 """return of string with (a subset of) the header lines 1651 (no terminal newline) 1652 """ 1653 # if no data, do not return a header 1654 if self.nvec == 0 or self.nt == 0: return '' 1655 1656 hlist = [] 1657 hlist.append('# ni_dimen = "%s"' % self.nt) 1658 if len(self.labels) == self.nvec: 1659 hlist.append('# ColumnLabels = "%s"' % ' ; '.join(self.labels)) 1660 if len(self.groups) == self.nvec: 1661 hlist.append('# ColumnGroups = "%s"' 1662 % self.val_list_to_at_str(self.groups)) 1663 hlist.append('# RowTR = "%s"' % self.tr) 1664 if len(self.goodlist) > 0: 1665 hlist.append('# GoodList = "%s"' % UTIL.encode_1D_ints(self.goodlist)) 1666 if self.nrowfull > 0: 1667 hlist.append('# NRowFull = "%s"' % self.nrowfull) 1668 if len(self.runstart) > 0: 1669 hlist.append('# RunStart = "%s"' \ 1670 % UTIL.int_list_string(self.runstart, sepstr=',')) 1671 1672 return '\n'.join(hlist) 1673 1674 def val_list_to_at_str(self, vlist): 1675 """convert a list like 2, 2, 2, 5, 9, 9, 11 to a string of the form 1676 3@2,5,2@9,11 1677 do not worry about types 1678 """ 1679 # start by making a vcount array 1680 vlen = len(vlist) 1681 vstart = 0 1682 slist = [] 1683 while vstart < vlen: 1684 val = vlist[vstart] 1685 vcur = vstart + 1 1686 while vcur < vlen and val == vlist[vcur]: 1687 vcur += 1 1688 nval = vcur-vstart 1689 if nval == 1: slist.append('%s' % val) 1690 else: slist.append('%d@%s' % (nval, val)) 1691 vstart = vcur 1692 return ','.join(slist) 1693 1694 def split_into_padded_runs(self, prefix, overwrite=0): 1695 """write one 1D file per run, where each is the same as the input, but 1696 is zero for every non-current run (i.e. only 1 non-zero run each) 1697 """ 1698 1699 if self.verb > 1: 1700 print('++ splitting into %d padded runs (rows=%d, cols=%d, prefix=%s)'\ 1701 % (self.nruns, self.nt, self.nvec, prefix)) 1702 1703 # if we don't have multiple runs, there is little to do 1704 if self.nruns <= 1: 1705 return self.write('%s.r01.1D' % prefix, overwrite=overwrite) 1706 1707 # store the data (self.mat) separately and start with 0-filled matrix 1708 orig_mat = self.mat 1709 self.mat = [[0 for row in range(self.nt)] for col in range(self.nvec)] 1710 1711 # per run: copy orig data, write, zero copied data data (back to all-0) 1712 end = 0 1713 for rind, rlen in enumerate(self.run_len): 1714 start = end 1715 end += self.run_len[rind] 1716 fname = '%s.r%02d.1D' % (prefix,rind+1) 1717 if self.verb>2: print('-- write %s, rows %d..%d' % (fname,start,end-1)) 1718 # copy orig data 1719 for col in range(self.nvec): 1720 for row in range(start, end): 1721 self.mat[col][row] = orig_mat[col][row] 1722 # write 1723 if self.write(fname, overwrite=overwrite): return 1 1724 # re-zero copied data 1725 for col in range(self.nvec): 1726 for row in range(start, end): 1727 self.mat[col][row] = 0 1728 1729 # and re-insert matrix 1730 del(self.mat) 1731 self.mat = orig_mat 1732 1733 return 0 1734 1735 def write_as_timing(self, fname, invert=0, column=0): 1736 """write set TRs to a timing file, one run per run, using '*' for 1737 empty runs (two for first run) 1738 1739 If multi-column data, can choose one. 1740 1741 fname : required filename 1742 invert : write times for zero data TRs 1743 column : which column to write times as 1744 1745 return status""" 1746 1747 if self.verb > 3: print('-- Afni1D write_timing to %s, inv=%d, col=%d' \ 1748 %(fname, invert, column)) 1749 1750 if not self.ready: 1751 print("** Afni1D not ready for write_timing to '%s'" % fname) 1752 return 1 1753 1754 err, tstr = UTIL.make_timing_string(self.mat[column], 1755 self.nruns, self.tr, invert) 1756 if err: return 1 1757 1758 if fname == '-' or fname == 'stdout': 1759 fp = sys.stdout 1760 else: 1761 try: fp = open(fname, 'r') 1762 except: 1763 print("** failed to open file '%s'" % fname) 1764 err = 1 1765 1766 if err: return 1 1767 1768 fp.write(tstr) 1769 if fp != sys.stdout: 1770 fp.close() 1771 1772 return 0 1773 1774 def write_censortr(self, fname, invert=0, column=0): 1775 """write set TRs to a timing file, one run per run, using '*' for 1776 empty runs (two for first run) 1777 1778 If multi-column data, can choose one. 1779 1780 fname : required filename 1781 invert : write times for zero data TRs 1782 column : which column to write times as 1783 1784 return status""" 1785 1786 if self.verb > 3: print('-- Afni1D write_censortr to %s, inv=%d, col=%d' \ 1787 %(fname, invert, column)) 1788 1789 if not self.ready: 1790 print("** Afni1D not ready for write_censortr to '%s'" % fname) 1791 return 1 1792 1793 if not self.run_info_is_consistent(whine=1): return 1 1794 1795 err,cstr = UTIL.make_CENSORTR_string(self.mat[column], 1796 rlens=self.run_len, asopt=1, verb=self.verb) 1797 1798 if err: return 1 1799 1800 try: fp = open(fname, 'w') 1801 except: 1802 print("** failed to open file '%s'" % fname) 1803 err = 1 1804 1805 if err: return 1 1806 1807 fp.write(cstr) 1808 fp.write('\n') # add a newline 1809 fp.close() 1810 1811 def append_vecs(self, matlist, newname=''): 1812 """append each Afni1D to the current one""" 1813 # test before trashing current matrix 1814 if not self.ready: 1815 print('** append: Afni1D is not ready') 1816 return 1 1817 1818 if self.verb > 3: print('-- Afni1D append_vecs...') 1819 1820 # allow matlist to be a simple mat 1821 if type(matlist) == type(self): 1822 newmats = [matlist] 1823 elif type(matlist) == type([]): 1824 if type(matlist[0]) == type(self): newmats = matlist 1825 else: 1826 print('** append_vecs: matlist elements not Afni1D') 1827 return 1 1828 else: 1829 print('** append_vecs: matlist must be list of Afni1D') 1830 return 1 1831 1832 for mat in newmats: 1833 if not mat.ready: 1834 print('** append: Afni1D is not ready') 1835 return 1 1836 if mat.nt != self.nt: 1837 print('** append: nt differs (%d != !%d)' % (mat.nt, self.nt)) 1838 return 2 1839 1840 # update labels 1841 empty = 1 1842 if not self.labels: self.labels = ['' for ind in range(self.nvec)] 1843 else: empty = 0 1844 for mat in newmats: 1845 if not mat.labels: 1846 self.labels.extend(['' for ind in range(mat.nvec)]) 1847 else: 1848 self.labels.extend(mat.labels) 1849 empty = 0 1850 if empty: 1851 del(self.labels) 1852 self.labels = [] 1853 1854 # actual appending... 1855 for newdset in newmats: 1856 self.mat.extend(newdset.mat) 1857 self.nvec += newdset.nvec 1858 1859 # update filename 1860 if newname: self.fname = newname 1861 elif self.fname.find('appended') < 0: self.fname += ' (appended)' 1862 1863 # nuke things that no longer apply 1864 del(self.groups) 1865 self.groups = [] 1866 del(self.goodlist) 1867 self.goodlist = [] 1868 del(self.basisinfo) 1869 self.basisinfo = [] 1870 del(self.cormat) 1871 del(self.cosmat) 1872 self.cormat = None 1873 self.cosmat = None 1874 self.cormat_ready = 0 1875 1876 return 0 1877 1878 def simplecopy(self): 1879 return Afni1D(from_mat=1, matrix=self.mat, verb=self.verb) 1880 1881 def show_group_labels(self): 1882 show_groups = (len(self.groups) == self.nvec) 1883 show_labs = (len(self.labels) == self.nvec) 1884 if not show_groups and not show_labs: 1885 print('** no label info to show') 1886 return 1887 1888 for ind in range(self.nvec): 1889 if self.verb: 1890 if show_groups: gstr = ', group %-3s' % self.groups[ind] 1891 else: gstr = '' 1892 if show_labs: lstr = ', label %s' % self.labels[ind] 1893 else: lstr = '' 1894 print('index %3d%s%s' % (ind, gstr, lstr)) 1895 elif show_labs: print('%s' % self.labels[ind]) 1896 1897 def show_labels(self): 1898 if self.verb: 1899 print('++ labels are:', self.labels) 1900 else: 1901 print(' '.join(self.labels)) 1902 1903 def show_major_order_of_labels(self): 1904 """be picky and verify that labels look like sSLICE.NAME, where 1905 SLICE increments (slowly) over slice index""" 1906 1907 if self.verb > 3: print('-- Afni1D show_major_order_of_labels...') 1908 1909 if not self.labels: 1910 print('** no labels to test for ordering') 1911 return 1912 1913 try: 1914 # nuke '.' and leading 's', remainder should be sorted ints 1915 vallist = [l.split('.')[0] for l in self.labels] 1916 vallist = [l.split('s')[1] for l in vallist] 1917 vallist = [int(i) for i in vallist] 1918 sorted = 1 1919 for ind in range(len(vallist)-1): 1920 if vallist[ind] > vallist[ind+1]: 1921 sorted = 0 1922 break 1923 if sorted: print('++ labels in row-major order: YES') 1924 else: print('++ labels in row-major order: NO') 1925 except: 1926 print('++ labels are not of expected slice-based format, '\ 1927 'cannot determine ordering') 1928 self.show_labels() 1929 1930 def clear_cormat(self): 1931 """nuke any existing cormat""" 1932 if self.verb > 3: print('-- Afni1D clear_cormat...') 1933 if not self.ready: return 1934 if self.cormat_ready: 1935 del(self.cormat) 1936 del(self.cosmat) 1937 self.cormat = None 1938 self.cosmat = None 1939 self.cormat_ready = 0 1940 1941 def set_cormat(self, update=0): 1942 """set cormat (the correlation matrix) and cormat.ready 1943 set cosmat (cosine matrix) as a follower 1944 1945 Note that the (Pearson's) correlations are de-meaned, and so the 1946 constant baseline terms can be highly correlated. 1947 1948 The cosine matrix might also be quite helpful. 1949 """ 1950 if self.verb > 3: print('-- Afni1D set_cormat...') 1951 if not self.ready: return 1952 if self.cormat_ready: 1953 if not update: return 1954 self.clear_cormat() # otherwise, nuke the old stuff and re-generate 1955 1956 if self.nvec < 2 or self.nt < 2: return 1957 1958 # make copy to abuse 1959 try: cmat = copy.deepcopy(self) 1960 except: 1961 print('... deepcopy failure, using simplecopy()...') 1962 cmat = self.simplecopy() 1963 1964 # demean each vector (for cormat), unless it is constant 1965 means = [UTIL.loc_sum(vec)/float(cmat.nt) for vec in cmat.mat] 1966 for v in range(cmat.nvec): 1967 lmin = min(cmat.mat[v]) 1968 lmax = max(cmat.mat[v]) 1969 # rcr - why avoid this and leave constant terms? 1970 if lmin != lmax: 1971 for ind in range(cmat.nt): 1972 cmat.mat[v][ind] -= means[v] 1973 1974 # and normalize 1975 norms = [UTIL.euclidean_norm(row) for row in cmat.mat] 1976 for v in range(cmat.nvec): 1977 for ind in range(cmat.nt): 1978 if norms[v] == 0.0: cmat.mat[v][ind] = 0.0 1979 else: cmat.mat[v][ind] /= norms[v] 1980 1981 # finally, assign cormat 1982 self.cormat =[[UTIL.dotprod(r1,r2) for r2 in cmat.mat] for r1 in cmat.mat] 1983 self.cormat_ready = 1 1984 1985 # and generate cosmat (dot product and scale) 1986 cnorm = [ UTIL.euclidean_norm(v) for v in self.mat ] 1987 cosmat = [[UTIL.dotprod(r1,r2) for r2 in self.mat] for r1 in self.mat] 1988 for v1 in range(len(cosmat)): 1989 for v2 in range(len(cosmat)): 1990 prod = cnorm[v1]*cnorm[v2] 1991 if prod != 0.0: cosmat[v1][v2] /= prod 1992 self.cosmat = cosmat 1993 1994 # nuke temporary (normalized) matrices 1995 del(cmat) 1996 del(means) 1997 del(norms) 1998 del(cnorm) 1999 2000 def show_cormat_diff_wout_gmean(self, unit=0, dp=3, spaces=2): 2001 ccopy = self.get_cormat_diff_wout_gmean(unit=unit, dp=dp, spaces=spaces) 2002 ccopy.show_cormat(dp=dp, spaces=spaces) 2003 del(ccopy) 2004 2005 def get_cormat_diff_wout_gmean(self, unit=0, dp=3, spaces=2): 2006 adcopy = self.copy() 2007 adcopy.project_out_vec(unit=unit) 2008 self.set_cormat(update=1) 2009 adcopy.set_cormat(update=1) 2010 2011 for rind, row in enumerate(adcopy.cormat): 2012 for ind in range(len(row)): 2013 row[ind] = self.cormat[rind][ind] - row[ind] 2014 return adcopy 2015 2016 def show_cormat(self, dp=3, spaces=2): 2017 """print the correlation matrix, to the given number of places 2018 dp > 0 : use that to right of decimal 2019 dp = 0 : use %g 2020 dp = -1 : use %f 2021 """ 2022 self.set_cormat() 2023 self.show_mat(self.cormat, dp=dp, spaces=spaces) 2024 2025 def show_distmat(self, dp=4, spaces=2, verb=None): 2026 """print a distance matrix, to the given number of places 2027 2028 Treating the input as lines of vectors (coordinates), print out an 2029 nrows x nrows distance matrix, the Euclidian distance between the 2030 pairs of coordinates. 2031 2032 verb: integral verbosity level (None : use self.verb) 2033 2034 return 0 on success 2035 """ 2036 if verb is None: 2037 verb = self.verb 2038 2039 if not self.ready: 2040 if verb: print('** no matrix to compute distance matrix from') 2041 return 1 2042 2043 if verb > 2: print("distmat: nvec = %d, nt = %d" % (self.nvec, self.nt)) 2044 2045 # get a transposed matrix, to assume coords are per row 2046 trdata = self 2047 trdata.transpose() 2048 vlen = trdata.nt 2049 2050 dmat = [] 2051 for cind, cvect in enumerate(trdata.mat): 2052 drow = [] 2053 for dind, dvect in enumerate(trdata.mat): 2054 diff = [cvect[i]-dvect[i] for i in range(vlen)] 2055 drow.append(UTIL.euclidean_norm(diff)) 2056 dmat.append(drow) 2057 2058 trdata.show_mat(dmat, dp=dp, spaces=spaces) 2059 2060 del(trdata) 2061 2062 return 0 2063 2064 def show_mat(self, mat, dp=3, spaces=2): 2065 """print the correlation matrix, to the given number of places 2066 dp > 0 : use that to right of decimal 2067 dp = 0 : use %g 2068 dp = -1 : use %f 2069 """ 2070 ss = ' '*spaces 2071 # add space for alignment to include a sign,0,. 2072 dt = dp + 3 2073 for v, vec in enumerate(mat): 2074 for val in vec: 2075 if dp > 0: ps = "%*.*f%s" % (dt, dp, val, ss) 2076 elif dp == 0: ps = "%g%s" % (val, ss) 2077 else: ps = "%f%s" % (val, ss) 2078 print(ps, end='') 2079 print("") 2080 2081 def make_cormat_warnings_string(self, cutoff=0.4, name='', 2082 skip_expected=1): 2083 """make a string for any entires at or above cutoffs: 2084 cut0=1.0, cut1=(1.0+cutoff)/2.0, cut2=cutoff 2085 2086 cut0, cut1, cut2 are cutoff levels (cut0=highest) 2087 that determine the severity of the correlation 2088 (anything below cut2 is ignored) 2089 2090 return error code (0=success) and 'warnings' string""" 2091 2092 if self.verb > 3: print("-- make_cormat_warn_str for '%s', cut=%g" \ 2093 % (name, cutoff)) 2094 2095 # assign base cutoff 2096 if cutoff < 0.0 or cutoff >= 1.0 : cutoff = 0.4 2097 2098 cut0 = 1.0 2099 cut1 = (1.0 + cutoff)/2.0 2100 cut2 = cutoff 2101 2102 # badlist holds cov, cos, row, col 2103 err, errstr, badlist = self.list_cormat_warnings(cutoff=cut2, 2104 skip_expected=skip_expected) 2105 if err: return err, errstr 2106 2107 blen = len(badlist) 2108 2109 if blen == 0: 2110 return 0, '-- no warnings for correlation matrix (cut = %.3f) --'%cut2 2111 2112 mstr = '\nWarnings regarding Correlation Matrix: %s\n\n' % name 2113 2114 # check for any constant 0 vectors here, too 2115 nfound = 0 2116 for ind in range(self.nvec): 2117 if UTIL.vals_are_constant(self.mat[ind], 0): 2118 mstr += ' ** vector #%02d is all ZEROs\n' % ind 2119 nfound += 1 2120 if nfound > 0: mstr += '\n' 2121 2122 mstr = mstr + \ 2123 ' severity correlation cosine regressor pair\n' \ 2124 ' -------- ----------- ------ ' + 40*'-' + '\n' 2125 2126 cutstrs = [ ' IDENTICAL: ', ' high: ', ' medium: ' ] 2127 2128 # note the maximum label length 2129 if self.labels: 2130 mlab = max([len(self.labels[col]) for val, s, row, col in badlist]) 2131 2132 for val, s, row, col in badlist: 2133 if abs(val) >= cut0: cs = cutstrs[0] 2134 elif abs(val) >= cut1: cs = cutstrs[1] 2135 else: cs = cutstrs[2] 2136 2137 # we have an appropriately evil entry... 2138 if self.labels: 2139 mstr += '%s %6.3f %6.3f (%2d vs. %2d) %*s vs. %s\n' % \ 2140 (cs, val, s, col, row, 2141 mlab, self.labels[col], self.labels[row]) 2142 else: 2143 mstr += '%s %6.3f %6.3f #%2d vs. #%2d\n' % \ 2144 (cs, val, s, col, row) 2145 2146 return 0, mstr 2147 2148 def list_cormat_warnings(self, cutoff=0.4, skip_expected=1): 2149 """return an error code, error string and a list of corval, cosval, 2150 vec, index for each cormat value with abs() > cutoff 2151 2152 if skip_expected, skip: mot or base against mot or base 2153 """ 2154 2155 if self.verb > 3: print('-- Afni1D list_cormat_warnings, cut=%g'%cutoff) 2156 2157 if not self.ready: 2158 return 1, '** no X-matrix to compute correlation matrix from', None 2159 2160 if not self.cormat_ready: self.set_cormat() # create cormat 2161 if not self.cormat_ready: # then failure 2162 return 1, '** cormat_warnings: failed to create cormat', None 2163 2164 cmat = self.cormat 2165 smat = self.cosmat 2166 2167 basecols = self.cols_by_group_list([-1]) 2168 motcols = self.cols_by_group_list([0]) 2169 roicols = self.cols_by_group_list([],allroi=1) 2170 2171 if self.verb > 1: 2172 print('-- LCP: len(base, mot, roi) = (%d, %d, %d), cut = %.2f' % \ 2173 (len(basecols), len(motcols), len(roicols), cutoff)) 2174 2175 # make a list of (abs(val),val,cosine,r,c) tuples in lower triangle 2176 clist = [] 2177 for r in range(1,self.nvec): 2178 for c in range(r): 2179 clist.append((abs(cmat[r][c]), cmat[r][c], smat[r][c], r, c)) 2180 2181 # clist.sort(reverse=True) # fails on old versions 2182 clist.sort() # smallest to largest, so process from end 2183 2184 # now make a list of evil-doers 2185 badlist = [] 2186 2187 # process list as smallest to largest, since old sort had no reverse 2188 clen = len(clist) 2189 for ind in range(clen): 2190 aval, val, s, r, c = clist[clen-1-ind] 2191 2192 # now raised to apqc, so do not auto-flag 2-run polort 0 correlations 2193 # [28 Oct 2021 rickr - noted by PT] 2194 # if aval == 1.0: 2195 # badlist.append((val, s, r, c)) # flag duplication 2196 # continue 2197 2198 # skip motion against either motion or baseline 2199 rbase = r in basecols 2200 cbase = c in basecols 2201 rmot = r in motcols 2202 cmot = c in motcols 2203 2204 # also, skip baseline against baseline (from 2-3 run polort 0) 2205 # 30 Dec 2013 2206 if skip_expected and (cmot or cbase) and (rbase or rmot): continue 2207 2208 if aval < cutoff: break 2209 2210 badlist.append((val, s, r, c)) # so keep this one 2211 2212 if self.verb > 1: 2213 print('-- LCP: badlist length = %d' % len(badlist)) 2214 2215 del(clist) 2216 2217 return 0, '', badlist 2218 2219 def labs_matching_str(self, mstr): 2220 if type(self.labels) != type([]): return [] 2221 return [lab for lab in self.labels if lab.find(mstr) >= 0] 2222 2223 def update_group_info(self): 2224 """if self.groups, try to set nruns and nroi""" 2225 2226 self.run_len = [self.nt] # init to 1 run 2227 2228 self.nroi = self.nvec # init to all 2229 if self.groups: 2230 if self.nvec: self.set_nruns() 2231 if len(self.groups) > 0: 2232 self.nroi = max(self.groups) 2233 if self.nroi < 0: self.nroi = 0 2234 2235 def set_nruns(self, nruns=0, run_lens=[]): 2236 """if nruns is positive, apply 2237 else if len(run_lens) > 0, apply 2238 else if have RunStart list, apply 2239 else, try to apply from a single run: 2240 find the first column in group 0, verify that it looks like 2241 one run of 1s (with remaining 0s), and use it to set the length 2242 2243 --> set nruns, run_len and runstart 2244 2245 return 0 on success, 1 on error 2246 """ 2247 2248 if self.verb > 3: 2249 print('-- LAD:set_nruns (nruns = %d, run_lens = %s)' % (nruns,run_lens)) 2250 print(' len(goodlist) = %d, runstart = %s' \ 2251 % (len(self.goodlist), self.runstart)) 2252 2253 # try to apply any passed nruns, first 2254 if nruns > 0: 2255 rlen = self.nt // nruns 2256 if rlen * nruns == self.nt: 2257 self.nruns = nruns 2258 self.run_len = [rlen for run in range(nruns)] 2259 2260 # update non-censored version 2261 if self.run_len_nc[0] == 0: 2262 self.run_len_nc = [0 for ll in self.run_len] 2263 elif self.run_len_nc[0] < self.nt: 2264 print('** cannot reset nruns in face of censored TRs') 2265 return 1 2266 else: 2267 self.run_len_nc = self.run_len[:] 2268 2269 if self.verb > 1: print('++ set_nruns: nruns = %d' % nruns) 2270 else: 2271 print('** nvalid nruns = %d (does not divide nt = %d)' \ 2272 % (nruns, self.nt)) 2273 return 1 2274 # and set runstart 2275 self.runstart = [r*rlen for r in range(nruns)] 2276 if self.verb > 1: 2277 print('++ set_nruns 0: nruns = %d, run_len = %s, runstart = %s' \ 2278 % (nruns, self.run_len, self.runstart)) 2279 return 0 2280 2281 # next, try run_lens (if not set, try runstart data from RunStart label) 2282 if type(run_lens) != type([]): 2283 print("** set_runs: run_lens is not a list type: %s" % run_lens) 2284 return 1 2285 2286 # if no run_lens but have self.runstart, convert and apply same logic 2287 if len(run_lens) == 0 and len(self.runstart) > 0: 2288 # first try to set run lengths taking goodlist into account 2289 if not self.apply_goodlist(): return 0 # all is well 2290 else: 2291 rlens = self.runstart[:] 2292 for ind in range(len(rlens)-1): 2293 rlens[ind] = rlens[ind+1] - rlens[ind] 2294 rlens[-1] = self.nt - rlens[-1] 2295 if self.verb > 1: 2296 print('++ setting run_len based on RunStart list:' , rlens) 2297 2298 rlens = run_lens 2299 2300 if len(rlens) > 0: 2301 if not UTIL.vals_are_positive(rlens): 2302 print("** set_runs: non-positive run length in list: %s" % rlens) 2303 return 1 2304 if UTIL.loc_sum(rlens) != self.nt: 2305 print("** set_runs: sum of run lengths != nt (%d): %s" \ 2306 % (self.nt, rlens)) 2307 return 1 2308 self.run_len = rlens[:] # make a copy, to be safe 2309 self.nruns = len(rlens) 2310 rsum = 0 2311 self.runstart = [] 2312 for ll in rlens: 2313 self.runstart.append(rsum) 2314 rsum += ll 2315 if self.verb > 1: 2316 print('++ set_nruns 1: nruns = %d, run_len = %s, runstart = %s' \ 2317 % (nruns, rlens, self.runstart)) 2318 return 0 2319 2320 # otherwise, init to 1 run and look for other labels 2321 self.nruns = 1 2322 if not self.groups or not self.nvec: return 0 2323 2324 # rcr : update poor use of baselines? 2325 errs = 0 2326 try: 2327 base_ind = self.groups.index(-1) 2328 if base_ind < 0: 2329 if self.verb > 1: print('-- no baseline group to set runs with') 2330 return 1 2331 2332 b0 = self.mat[base_ind] 2333 2334 # skip 0s, count 1s, rest must be 0s 2335 for val in b0: 2336 if val != 0 and val != 1: 2337 if self.verb > 1: print('-- baseline vals not just 0,1: %s' % val) 2338 return 1 2339 2340 # looks good, find run length and the number of runs 2341 first = b0.index(1) # find first 1 2342 try: next = b0.index(0,first+1) # find next 0 2343 except: next = self.nt 2344 if next <= first: 2345 if self.verb > 1: print('-- odd run_len check...') 2346 return 1 2347 2348 # we have a run length 2349 rlen = next - first 2350 if rlen > 0: nruns = self.nt//rlen # integral division 2351 else: 2352 if self.verb > 1: 2353 print('** failed to set rlen from baseline (%d,%d)'%(first,next)) 2354 return 1 2355 if rlen*nruns != self.nt: 2356 if self.verb>1: 2357 print('** nruns failure: rlen %d, nruns %d, len %d' % \ 2358 (rlen, nruns, len(b0))) 2359 if self.nrowfull > self.nt: 2360 print('++ nrowfull (%d) > nt (%d) ==> censored TRs' \ 2361 % (self.nrowfull, self.nt)) 2362 return 1 2363 2364 # success! 2365 2366 self.run_len = [rlen for i in range(nruns)] 2367 self.nruns = nruns 2368 2369 except: 2370 if self.verb > 1: print('** unknown exception in LD:set_nruns') 2371 errs = 1 2372 2373 return errs 2374 2375 def apply_goodlist(self, parent=None, padbad=0): 2376 """try to apply goodlist, runstart, nt and nrowfull 2377 -> set nruns and run_lens[] 2378 -> also, set run_lens_nc[] (not censored) 2379 2380 if padbad is set, pad the data with zero over all TRs NOT 2381 in goodlist (so run_lens should come from RunStart) 2382 2383 if parent is set (as an Afni1D instance), 2384 use it for goodlist, runstart, nrowfull 2385 2386 return 0 on success""" 2387 2388 if isinstance(parent, Afni1D): 2389 if self.verb > 2: print('-- apply_goodlist: run info from parent') 2390 2391 # allow for processing simple 1D 2392 if not parent.havelabs: 2393 return self.uncensor_from_vector(parent.mat[0]) 2394 2395 runstart = parent.runstart 2396 goodlist = parent.goodlist 2397 nrowfull = parent.nrowfull 2398 else: 2399 if self.verb > 2: print('-- apply_goodlist: run info from self') 2400 runstart = self.runstart 2401 goodlist = self.goodlist 2402 nrowfull = self.nrowfull 2403 2404 # first see if we have the necessary fields 2405 rv = 1 2406 try: 2407 if len(runstart) > 0 and len(goodlist) > 0 and \ 2408 self.nt > 0 and nrowfull > 0: rv = 0 2409 except: 2410 if self.verb > 1: print('** bad internals for apply_goodlist') 2411 2412 if rv == 1: return 1 # bail 2413 2414 # other tests 2415 if not UTIL.vals_are_sorted(goodlist): 2416 if self.verb > 1: print('** LAD: goodlist not sorted') 2417 return 1 2418 if len(goodlist) != self.nt: 2419 if self.verb > 1: print('** LAD: goodlist length (%d) != nt (%d)' \ 2420 % (len(goodlist), self.nt)) 2421 return 1 2422 if not UTIL.vals_are_sorted(runstart): 2423 if self.verb > 1: print('** LAD: runstart not sorted') 2424 return 1 2425 if goodlist[-1] >= nrowfull: 2426 if self.verb > 1: print('** LAD: max goodlist value exceeds rowfull') 2427 return 1 2428 if goodlist[0] < 0: 2429 if self.verb > 1: print('** LAD: goodlist has negative value') 2430 return 1 2431 if runstart[-1] >= nrowfull: 2432 if self.verb > 1: print('** LAD: max runstart value exceeds rowfull') 2433 return 1 2434 if runstart[0] < 0: 2435 if self.verb > 1: print('** LAD: runstart has negative value') 2436 return 1 2437 2438 # ---- now try to sort things out ---- 2439 2440 # if padding 'bad' TRs run_len is basic, but mat needs to be zero-filled 2441 if padbad: 2442 if self.GLapplied == 2: 2443 if self.verb > 0: print('** padbad: run already padded!') 2444 return 1 2445 2446 # first note the run lengths 2447 rlens = runstart[:] 2448 for ind in range(len(rlens)-1): 2449 rlens[ind] = rlens[ind+1] - rlens[ind] 2450 rlens[-1] = nrowfull - rlens[-1] 2451 2452 # then pad the data (create zero vectors, and copy good data) 2453 for vind in range(len(self.mat)): 2454 oldvec = self.mat[vind] 2455 newvec = [0 for t in range(nrowfull)] 2456 for gind in range(len(goodlist)): 2457 newvec[goodlist[gind]] = oldvec[gind] 2458 self.mat[vind] = newvec 2459 del(oldvec) 2460 2461 self.nt = nrowfull # now the lists are full 2462 self.nruns = len(runstart) # 2463 self.run_len = rlens # steal reference 2464 self.GLapplied = 2 # goodlist was applied 2465 2466 return 0 2467 2468 # --- not zero padding, so adjust run lengths --- 2469 2470 # set rnext list to be index *after* the current run 2471 rnext = runstart[1:] 2472 rnext.append(nrowfull) 2473 if self.verb > 3: print('-- set rnext list to %s' % rnext) 2474 nruns = len(runstart) 2475 2476 # accumulate run counts over goodlist 2477 run = 0 2478 rcount = [0 for r in runstart] # count entries per run 2479 for gval in goodlist: 2480 # maybe adjust the run index (note that rnext[-1] > goodlist vals) 2481 while gval >= rnext[run]: run += 1 2482 rcount[run] += 1 2483 2484 # and verify that we have accounted for everything 2485 gtot = UTIL.loc_sum(rcount) 2486 glen = len(goodlist) 2487 if gtot != glen: 2488 print('** apply_goodlist: gtot error: %d != %d' % (gtot, glen)) 2489 return 1 2490 if gtot != self.nt: 2491 print('** apply_goodlist: nt (%d) != goodtot (%d)' % (self.nt, gtot)) 2492 return 1 2493 2494 # also, generate an uncensored run length list from runstart 2495 self.run_len_nc = [runstart[i+1]-runstart[i] for i in range(nruns-1)] 2496 self.run_len_nc.append(self.nrowfull-self.runstart[-1]) 2497 2498 if self.verb > 1: 2499 print('++ from apply_goodlist: run_len[%d] = %s' % (nruns, rcount)) 2500 print(' run_len_nc = %s' % self.run_len_nc) 2501 print(' run_start = %s' % self.runstart) 2502 2503 self.nruns = nruns 2504 self.run_len = rcount # steal reference 2505 self.GLapplied = 1 # goodlist was applied 2506 2507 return 0 2508 2509 def show(self): 2510 print(self.make_show_str()) 2511 2512 def show_rows_cols(self, mesg='', verb=1): 2513 """display the number of rows (nt) and columns (nvec) 2514 2515 mesg: if set print before output 2516 verb: if 0, no text""" 2517 2518 if mesg: print('%s' % mesg, end='') 2519 if verb > 0: print('rows = %d, cols = %d' % (self.nt, self.nvec)) 2520 else: print('%d %d' % (self.nt, self.nvec)) 2521 2522 def get_tr_counts(self): 2523 """return status, self.run_len, self.run_len_nc 2524 2525 check for consistency 2526 if vector lengths are not nruns and 2527 """ 2528 2529 rv = self.check_tr_count_consistency(self.verb) 2530 return rv, self.run_len, self.run_len_nc 2531 2532 def check_tr_count_consistency(self, verb=1): 2533 """check the consistency of TR counts 2534 2535 - matrix is ready 2536 - nruns >= 1 2537 - len(run_len) == len(run_len_nc) == nruns 2538 - sum(run_len) == nt 2539 - sun(run_len_nc) == nrowfull 2540 2541 return 0 if consistent, 1 otherwise 2542 """ 2543 2544 if not self.ready: 2545 if verb: print('** bad run lists: not ready') 2546 return 1 2547 2548 nr = self.nruns 2549 if nr < 1: 2550 if verb: print('** bad run lists: nruns = %d' % nr) 2551 return 1 2552 2553 if len(self.run_len) != nr: 2554 if verb: print('** bad run lists: len(run_len) = %d, nruns = %d' \ 2555 % (len(self.run_len), self.nruns)) 2556 return 1 2557 2558 if len(self.run_len_nc) != nr: 2559 if verb: print('** bad run lists: len(run_len_nc) = %d, nruns = %d' \ 2560 % (len(self.run_len_nc), self.nruns)) 2561 return 1 2562 2563 ntr = UTIL.loc_sum(self.run_len) 2564 if ntr != self.nt: 2565 if verb: print('** bad run lists: sum(run_len) = %d, nt = %d' \ 2566 % (ntr, self.nt)) 2567 return 1 2568 2569 ntr = UTIL.loc_sum(self.run_len_nc) 2570 if ntr != self.nrowfull: 2571 if verb: print('** bad run lists: sum(run_len_nc) = %d, nrf = %d' \ 2572 % (ntr, self.nrowfull)) 2573 return 1 2574 2575 return 0 2576 2577 def apply_censor_dset(self, cset=None): 2578 """remove censored TRs from self.mat 2579 2580 ** for now, it becomes 1 combined run 2581 """ 2582 2583 mat = self.mat 2584 nt = self.nt 2585 2586 if self.verb > 1: print('-- censoring input data...') 2587 2588 if cset == None: return 2589 if cset.nvec == 0: return 2590 if self.nvec == 0: return 2591 2592 if cset.nt != nt: 2593 print('** ACD: censort length (%d) does not match NT (%d)' \ 2594 % (cset.nt, nt)) 2595 return 2596 2597 clist = cset.mat[0] # get censor list for ease of typing 2598 2599 if self.verb > 2: 2600 nkeep = UTIL.loc_sum([v for v in clist if v]) 2601 print('-- censoring from length %d to length %d' % (nt, nkeep)) 2602 2603 # if nothing is censored out, we're done 2604 if UTIL.vals_are_constant(clist, 1): return 2605 2606 # else, make a new matrix 2607 cmat = [[col[ind] for ind in range(nt) if clist[ind]] for col in mat] 2608 2609 del(self.mat) 2610 self.mat = cmat 2611 2612 # now clean up some extra fields 2613 2614 self.nt = len(cmat[0]) 2615 self.nrowfull = self.nt 2616 2617 # could deal with run lengths here, but save for later 2618 self.nruns = 1 2619 self.run_len = [self.nt] 2620 self.goodlist = [] 2621 self.runstart = [] 2622 2623 del(self.cormat) 2624 del(self.cosmat) 2625 self.cormat = [] 2626 self.cosmat = [] 2627 self.cormat_ready = 0 2628 2629 return 2630 2631 def get_censored_mat(self, cset=None): 2632 """return self.mat, restricted to uncensored TRs""" 2633 2634 mat = self.mat 2635 nt = self.nt 2636 2637 if cset == None: return mat 2638 if cset.nvec == 0: return mat 2639 2640 if cset.nt != nt: 2641 print('** GCM: censort length (%d) does not match NT (%d)' \ 2642 % (cset.nt, nt)) 2643 return None 2644 2645 clist = cset.mat[0] # get censor list for ease of typing 2646 2647 # if nothing is censored out, return the original 2648 if UTIL.vals_are_constant(clist, 1): return mat 2649 2650 # else, make a copy 2651 return [[col[ind] for ind in range(nt) if clist[ind]] for col in mat] 2652 2653 def get_max_displacement(self, dtype=1, cset=None): 2654 """compute the maximum pairwise displacement among the 2655 N-dimensional coordinates (e.g. motion params) 2656 2657 dtype = 1: enorm 2658 2659 - so compute max(dist(x,y)) over all x, y coordinate pairs 2660 2661 dtype = 0: max abs 2662 2663 compute max absolue pairwise displacement among the N dimensions 2664 2665 - so compute max(abs(v_i,j - v_i,k)) over all all dimensions i 2666 and j,k index pairs 2667 - or compute as max(max(v_i)-min(v_i)) over dimensions i 2668 2669 cset is an optional censor dataset (1=keep, 0=censor) 2670 """ 2671 2672 maxdiff = 0 2673 2674 mat = self.get_censored_mat(cset) 2675 if mat == None: mat = self.mat 2676 2677 if self.nvec > 0: nt = len(mat[0]) 2678 else: nt = 0 2679 2680 if dtype == 0: # max vector displacement 2681 maxdiff = max([max(vec)-min(vec) for vec in mat]) 2682 2683 else: # max euclidean distance 2684 2685 # get all distances, a 2D list of the form [[DIST, i, j]] 2686 # (so length is NT choose 2 = (NT*(NT-1)/2), which is possibly long) 2687 2688 en = UTIL.euclidean_norm 2689 dlist = [[en([mat[i][r1]-mat[i][r2] for i in range(self.nvec)]), 2690 r1, r2] for r1 in range(nt) for r2 in range(r1)] 2691 2692 # and grab the biggest 2693 dlist.sort() 2694 maxdiff = dlist[-1][0] 2695 2696 if self.verb > 1: 2697 i = dlist[-1][1] 2698 j = dlist[-1][2] 2699 print('++ max edisp %.10f between indices %d and %d' \ 2700 % (maxdiff, i, j)) 2701 2702 del(dlist) 2703 2704 return maxdiff 2705 2706 def get_max_displacement_str(self, mesg='', enorm=2, verb=1): 2707 """display the maximum displacement based on enorm 2708 enorm = 1 : absolute displacement over dims only 2709 enorm = 2 : enorm only 2710 enorm = 3 : both 2711 """ 2712 2713 rstr = '' # return string 2714 2715 if enorm == 1 or enorm == 3: 2716 maxabval = self.get_max_displacement(0) 2717 if enorm: 2718 maxenorm = self.get_max_displacement(1) 2719 2720 if verb > 0: 2721 if enorm == 1: rstr = 'max_param_diff = %g' % maxabval 2722 elif enorm == 3: rstr = 'max_param_diff = %g, max_enorm_diff = %g' \ 2723 % (maxabval, maxenorm) 2724 else: rstr = 'max_enorm_diff = %g' % maxenorm 2725 else: 2726 if enorm == 1: rstr = '%g' % maxabval 2727 elif enorm == 3: rstr = '%g %g' % (maxabval, maxenorm) 2728 else: rstr = '%g' % maxenorm 2729 2730 if mesg: return "%s %s" % (mesg, rstr) 2731 else: return rstr 2732 2733 def show_min_mean_max_stdev(self, col=-1, verb=1): 2734 """show min, mean, max, stdev for each column (unless col specified)""" 2735 2736 if verb: print("file %s (len %d)" % (self.fname, self.nt)) 2737 for cind, col in enumerate(self.mat): 2738 if verb: 2739 ps = " col %d: " % cind 2740 form = "min = %7.4f, mean = %7.4f, max = %7.4f, stdev = %7.4f" 2741 else: 2742 ps = '' 2743 form = "%7.4f %7.4f %7.4f %7.4f" 2744 print(ps + form % UTIL.min_mean_max_stdev(col)) 2745 2746 def show_trs_to_zero(self, col=-1, verb=1): 2747 """show number of TRs until constant zero 2748 (i.e. search from end for first non-zero value) 2749 """ 2750 2751 if verb: print("file %s (len %d)" % (self.fname, self.nt)) 2752 for cind, col in enumerate(self.mat): 2753 # set mind = index of last non-zero element, and increment 2754 mind = -1 2755 for ind in range(len(col)-1,-1,-1): 2756 if col[ind] != 0: 2757 mind = ind 2758 break 2759 mind += 1 # becomes a count to zero 2760 2761 if verb: print(" col %d: response length = %d" % (cind, mind)) 2762 else: print("%d" % mind, end='') 2763 print('') 2764 2765 def show_xmat_stim_info(self, label='', verb=1): 2766 """display Basis* field information from an xmat 2767 Name : stim label 2768 Option : -stim_* option name 2769 Formula : 3dDeconvolve basis function 2770 Columns : column selector (test, but do not convert) 2771 2772 if label: only show matching classes 2773 """ 2774 if label == '': 2775 print("-- have %d Basis elements:\n" % len(self.basisinfo)) 2776 2777 for bdict in self.basisinfo: 2778 keys = list(bdict.keys()) 2779 2780 # maybe we want to skip this (probably most) 2781 if label and label != 'ALL': 2782 if 'Name' not in keys: 2783 continue 2784 if bdict['Name'] != label: 2785 continue 2786 2787 if 'column_list' in keys: 2788 cols = bdict['column_list'] 2789 ncol = len(cols) 2790 else: 2791 ncol = 0 2792 for key in g_xmat_basis_labels: 2793 if key in keys: val = bdict[key] 2794 else: val = '' 2795 # if there is only 1 column, do not use .. notation 2796 if key == 'Columns' and ncol == 1: 2797 val = str(cols[0]) 2798 print("%-10s : %s" % (key, val)) 2799 print() 2800 2801 def get_xmat_stype_cols(self, stypes=['ALL']): 2802 """return columns for the given stim types, based on BasisOption entries 2803 (requires them to be stored in dict['column_list'] 2804 (i.e. cannot include baseline or -stim_file regressor columns) 2805 2806 times: -stim_times 2807 AM: -stim_times_AM1 or -stim_times_AM2 2808 AM1: -stim_times_AM1 2809 AM2: -stim_times_AM2 2810 IM: -stim_times_IM 2811 2812 """ 2813 # make a list of acceptible options 2814 if 'ALL' in stypes: 2815 optlist = ['-stim_times', '-stim_times_AM1', '-stim_times_AM2', 2816 '-stim_times_IM'] 2817 else: 2818 # dupes are not a problem, if the user happens to allow them 2819 optlist = [] 2820 if 'times' in stypes: 2821 optlist.extend(['-stim_times']) 2822 if 'AM' in stypes: 2823 optlist.extend(['-stim_times_AM1', '-stim_times_AM2']) 2824 if 'AM1' in stypes: 2825 optlist.extend(['-stim_times_AM1']) 2826 if 'AM2' in stypes: 2827 optlist.extend(['-stim_times_AM2']) 2828 if 'IM' in stypes: 2829 optlist.extend(['-stim_times_IM']) 2830 2831 if self.verb > 2: 2832 print("-- show_xmat_stype_cols convert: %s\n" \ 2833 " to: %s" % (', '.join(stypes), ', '.join(optlist))) 2834 2835 # for each Basis entry, if the Option is desirable, append the cols 2836 collist = [] 2837 for bdict in self.basisinfo: 2838 keys = list(bdict.keys()) 2839 2840 if 'Option' not in keys or 'column_list' not in keys: 2841 continue 2842 2843 sopt = bdict['Option'] 2844 scols = bdict['column_list'] 2845 2846 # does this have an option of interest? 2847 if sopt not in optlist: 2848 continue 2849 2850 if len(scols) == 0: 2851 print("** xmat_type_cols: no 'column_list' for stim opt %s"%sopt) 2852 return [] 2853 2854 # we have something, append do the list 2855 collist.extend(scols) 2856 2857 collist.sort() 2858 2859 return collist 2860 2861 def slice_order_to_times(self, verb=1): 2862 """given a slice order, resort index list to slice times 2863 2864 e.g. If TR=2 and the slice order is: 0 2 4 6 8 1 3 5 7 9 2865 Then the slices/times ordered by time (as input) are: 2866 2867 slices: 0 2 4 6 8 1 3 5 7 9 2868 times: 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2869 2870 And the slices/times ordered instead by slice index are: 2871 2872 slices: 0 1 2 3 4 5 6 7 8 9 2873 times: 0.0 1.0 0.2 1.2 0.4 1.4 0.6 1.6 0.8 1.8 2874 2875 It is this final list of times that is output. 2876 """ 2877 2878 # maybe we need to transpose 2879 trans = 0 2880 if self.nt == 1 and self.nvec > 1: 2881 self.transpose() 2882 trans = 1 2883 2884 # create a matrix of index/slice time pairs 2885 data = [[val,0] for val in self.mat[0]] 2886 for ind, val in enumerate(data): 2887 val[1] = self.tr * ind * 1.0 / self.nt 2888 2889 # sort on the indices 2890 data.sort() 2891 2892 # and return the times 2893 self.mat[0] = [val[1] for val in data] 2894 2895 if trans: self.transpose() 2896 2897 return 0 2898 2899 def add_offset(self, offset): 2900 """add offset value to every value in matrix""" 2901 2902 if not self.ready: 2903 print('** add_affset: data not ready') 2904 return 1 2905 2906 for row in self.mat: 2907 for ind in range(len(row)): row[ind] += offset 2908 2909 return 0 2910 2911 def rank(self, style='dense', reverse=0, base1=0, verb=1): 2912 """convert to the min to max rank 2913 i.e. for each value, apply its ordered index 2914 e.g. 3.4 -0.3 4.9 2.0 ==> 2 0 3 1 2915 2916 style : rank style ('dense' or 'competition') 2917 reverse: sort largest to smallest 2918 base1: apply 1-based ranks 2919 2920 If there is one row or col, use it. Otherwise, use column 0. 2921 2922 return status (0=success) 2923 """ 2924 2925 if not self.ready: 2926 print('** rank: data not ready') 2927 return 1, [] 2928 2929 if self.nvec == 0 or self.nt == 0: return 0 # nothing to do 2930 2931 # if nt == 1, use first value per vector (apply as transpose) 2932 if self.nt == 1 and self.nvec > 1: 2933 data = [self.mat[ind][0] for ind in range(self.nvec)] 2934 rv, data = UTIL.get_rank(data, style=style, reverse=reverse) 2935 if rv: return 1 2936 for ind in range(self.nvec): 2937 self.mat[ind][0] = data[ind] 2938 2939 # else nt > 1, process all vectors 2940 else: 2941 for ind, row in enumerate(self.mat): 2942 rv, data = UTIL.get_rank(row, style=style, reverse=reverse) 2943 if rv: return 1 2944 self.mat[ind] = data 2945 2946 if base1: return self.add_offset(1) 2947 else: return 0 2948 2949 def get_indices_str(self, ind_types): 2950 """return an index list (sub-brick selector form) for the 2951 following groups: 2952 2953 1: baseline (group -1) 2954 2: motion (group 0) 2955 4: interest (group > 0) 2956 2957 Do not return an empty list. If the groups do not exist or are 2958 not found, return '0..$'.""" 2959 2960 default = '0..$' 2961 2962 if self.verb > 1: print('-- show indices, types = %d, groups = %s' \ 2963 % (ind_types, self.groups)) 2964 2965 bmask = ind_types & 7 2966 if not self.ready: return default 2967 if bmask == 0 or bmask == 7: return default 2968 if len(self.groups) < 1: return default 2969 2970 ilist = [] 2971 allind = list(range(len(self.groups))) 2972 if ind_types & 1: 2973 ilist += [ind for ind in allind if self.groups[ind] == -1] 2974 if ind_types & 2: 2975 ilist += [ind for ind in allind if self.groups[ind] == 0] 2976 if ind_types & 4: 2977 ilist += [ind for ind in allind if self.groups[ind] > 0] 2978 ilist.sort() 2979 2980 elist = UTIL.encode_1D_ints(ilist) 2981 if elist == '': 2982 if self.verb > 1: print('-- replacing empty encode list with full') 2983 elist = default 2984 return elist 2985 2986 def make_show_str(self): 2987 if self.ready: rstr = 'ready' 2988 else: rstr = 'not ready' 2989 2990 mstr = "++ mat name : %s (%s)\n" \ 2991 "++ fname : %s\n" \ 2992 "++ nvec : %d\n" \ 2993 "++ nt : %d\n" \ 2994 "++ labels : %s\n" \ 2995 "++ groups : %s\n" \ 2996 "++ goodlist : %s\n" \ 2997 "++ runstart : %s\n" \ 2998 "++ tr : %s\n" \ 2999 "++ nrowfull : %d\n" \ 3000 "++ nruns : %d\n" \ 3001 "++ run_len : %s\n" \ 3002 "++ run_len_nc: %s\n" % \ 3003 (self.name, rstr, self.fname, self.nvec, self.nt, 3004 self.labels, self.groups, self.goodlist, self.runstart, 3005 self.tr, self.nrowfull, 3006 self.nruns, self.run_len, self.run_len_nc) 3007 3008 return mstr 3009 3010 def cols_by_group_list(self, groups, allroi=0): 3011 """return a list of columns, given a list of groups 3012 - if 'allroi' is 1, append to 'groups' all that are positive 3013 - if 'allroi' is 1, then 'group' must not contain any positive 3014 group indices""" 3015 if not self.groups: return [] 3016 3017 # append positive groups, but fail if some already exist 3018 if allroi: 3019 for val in groups: 3020 if val > 0: 3021 print("** CBGL: cannot pass positive groups with 'allroi'") 3022 return [] 3023 groups.extend([g for g in self.groups if g > 0]) 3024 if len(groups) < 1 or len(self.groups) < 1: return [] 3025 return [val for val in range(self.nvec) if self.groups[val] in groups] 3026 3027 def ordered_cols_by_group_list(self, groups): 3028 """return a list of columns, given a list of groups 3029 3030 for each group, insert all columns from that group 3031 (i.e. keep group order) 3032 """ 3033 if not self.groups: return [] 3034 3035 if len(groups) < 1 or len(self.groups) < 1: return [] 3036 3037 clist = [] 3038 for g in groups: 3039 clist.extend([v for v in range(self.nvec) if self.groups[v] == g]) 3040 3041 return clist 3042 3043 def cols_by_label_list(self, labels): 3044 """return a list of columns, given a list of labels""" 3045 if not self.labels or not labels: return [] 3046 if not list2_is_in_list1(self.labels, labels, "labels"): return [] 3047 # make a working function, suggested by I Schwabacher 3048 return [self.labels.index(lab) for lab in labels] 3049 3050 def init_from_matrix(self, matrix): 3051 """initialize Afni1D from a 2D (or 1D) array""" 3052 3053 if self.verb > 3: print("-- Afni1D: init_from_matrix") 3054 3055 if type(matrix) != type([]): 3056 print('** matrix for init must be [[float]] format') 3057 return 1 3058 elif type(matrix[0]) != type([]): 3059 if type(matrix[0]) == type(1.0): 3060 3061 # init from vector 3062 self.mat = [matrix] 3063 self.nvec = 1 3064 self.nt = len(matrix) 3065 self.ready = 1 3066 return 1 3067 3068 else: 3069 print('** matrix for init must be [[float]] format') 3070 return 1 3071 3072 self.mat = matrix 3073 self.nvec = len(matrix) 3074 self.nt = len(matrix[0]) 3075 self.ready = 1 3076 3077 return 0 3078 3079 def init_from_general_name(self, name): 3080 """might contain [] or {} selection 3081 might be '-' or stdin""" 3082 3083 aname = BASE.afni_name(self.fname) 3084 3085 if self.verb > 3: print("-- Afni1D: init_from_gen name '%s'" % name) 3086 3087 if self.init_from_1D(aname.rpve()): return 1 # failure 3088 3089 self.aname = aname 3090 self.fname = aname.rpve() 3091 self.name = aname.pve() 3092 3093 # apply column and/or row selectors 3094 if aname.colsel: 3095 if self.nvec == len(self.labels): 3096 labels = self.labels 3097 else: 3098 labels = [] 3099 3100 ilist = UTIL.decode_1D_ints(aname.colsel, verb=self.verb, 3101 imax=self.nvec-1, labels=labels) 3102 if ilist == None: return 1 3103 if self.reduce_by_vec_list(ilist): return 1 3104 if aname.rowsel: 3105 ilist = UTIL.decode_1D_ints(aname.rowsel,verb=self.verb,imax=self.nt-1) 3106 if ilist == None: return 1 3107 if self.reduce_by_tlist(ilist): return 1 3108 3109 return 0 3110 3111 def add_basisinfo(self, label, data): 3112 """update self.basisinfo 3113 handle basis types: 3114 Name : stim label 3115 Option : -stim_* option name 3116 Formula : 3dDeconvolve basis function 3117 Columns : column selector (test, but do not convert) 3118 3119 column_list : included for valid Columns 3120 """ 3121 if self.verb > 3: print("-- add_basisinfo: %s = %s" %(label,data)) 3122 3123 status, btype, bind = self.basisinfo_label_parts(label) 3124 if status < 0: 3125 print("** illegal basisinfo: %s = '%s'" % (label, data)) 3126 return 1 3127 elif status == 0: 3128 return 0 3129 3130 # we have something to process 3131 3132 # make sure there are enough dictionaries, and note the current one 3133 nprev = len(self.basisinfo) 3134 if nprev < bind: 3135 self.basisinfo.extend([{} for i in range(bind-nprev)]) 3136 bdict = self.basisinfo[bind-1] 3137 3138 # test data for given btype 3139 if btype == 'Name': 3140 pass 3141 elif btype == 'Option': 3142 if not data.startswith('-stim_'): 3143 if self.verb > 1: 3144 print("** basisinfo: bad btype in label for: %s = '%s'" \ 3145 % (label, data)) 3146 # what to do, just whine? 3147 # return 1 3148 elif btype == 'Formula': 3149 # consider setting things like stim dur, etc 3150 pass 3151 elif btype == 'Columns': 3152 # if we have a valid column list, store it 3153 # ** note: replace ':' with '..' 3154 dotdata = data.replace(':', '..') 3155 intlist = UTIL.decode_1D_ints(dotdata) 3156 if len(intlist) == 0: 3157 if self.verb > 1: 3158 print("** basisinfo: bad Columns in: %s = '%s'"%(label, data)) 3159 return 1 3160 3161 if self.nvec > 0: imax = self.nvec-1 3162 else: imax = -1 3163 bdict['column_list'] = intlist 3164 # store in .. format, not : 3165 data = dotdata 3166 3167 bdict[btype] = data 3168 # make a list of 3169 3170 return 0 3171 3172 def basisinfo_label_parts(self, label, prefix='Basis'): 3173 """return a tuple of: 3174 status, basistype, index 3175 3176 status = -1 (on error), 0 (on ignore), 1 (on success) 3177 """ 3178 lstuff = label.split('_') 3179 3180 # skip anything not of the format Basis*_index (1-based) 3181 if len(lstuff) < 2: return 0, '', 0 3182 3183 blen = len(prefix) 3184 btype = lstuff[0][blen:] 3185 bind = -1 3186 try: 3187 bind = int(lstuff[1]) 3188 except: 3189 print("** bad baseinfo label: %s") 3190 return -1, '', 0 3191 3192 return 1, btype, bind 3193 3194 def init_from_1D(self, fname): 3195 """initialize Afni1D from a 1D file (return err code)""" 3196 3197 if self.verb > 3: print("-- Afni1D: init_from_1D '%s'" % fname) 3198 3199 tmat, clines = TD.read_data_file(fname, verb=self.verb) 3200 if not TD.data_is_rect(tmat): 3201 print("** data is not rectangular in %s" % fname) 3202 return 1 3203 # if not tmat: return 1 3204 3205 # treat columns as across time 3206 mat = UTIL.transpose(tmat) 3207 del(tmat) 3208 3209 self.mat = mat 3210 self.nvec = len(mat) 3211 if self.nvec > 0: self.nt = len(mat[0]) 3212 else: self.nt = 0 3213 self.ready = 1 3214 3215 # last step, process comment lines 3216 if not clines: return 0 3217 3218 for line in clines: 3219 label, data = c1D_line2labelNdata(line) 3220 if not label: 3221 # store all unprocessed comment lines 16 Apr 2018 3222 self.header.append(line) 3223 continue 3224 3225 verb_level = 3 # cutoff for verbose level 3226 3227 self.havelabs = 1 3228 3229 try: # to process some of the data 3230 if label == 'ni_type': 3231 if data.find('*') >= 0: 3232 ncols, dtype = data.split('*') 3233 else: 3234 ncols = 1 3235 dtype = data 3236 ncols = int(ncols) 3237 if self.verb > verb_level: 3238 print("-- label %s: cols = %d, dtype = %s" % \ 3239 (label, ncols, dtype)) 3240 if ncols != self.nvec: 3241 print("** matrix vecs %d != %s cols %d" % \ 3242 (self.nvec, label, ncols)) 3243 elif label == 'ni_dimen': 3244 nrows = int(data.split(',')[0]) # 3D datasets have commas 3245 if self.verb > verb_level: 3246 print("-- label %s: rows = %d" % (label, nrows)) 3247 if nrows != self.nt: 3248 print("** matrix nt %d != %s rows %d" % \ 3249 (self.nt, label, nrows)) 3250 elif label == 'ColumnLabels': 3251 if data == '': 3252 self.labels = [] 3253 else: 3254 self.labels = [ss.strip() for ss in data.split(';')] 3255 if self.verb > verb_level: 3256 print("-- label %s: labels = %s" % (label, self.labels)) 3257 if self.nvec != len(self.labels): 3258 print("** %d ColumnLabels but %d columns" % \ 3259 (len(self.labels), self.nvec)) 3260 self.labels = [] 3261 elif label == 'ColumnGroups': 3262 self.groups = UTIL.decode_1D_ints(data) 3263 if self.groups: 3264 if len(self.groups) != self.nvec: 3265 print("** ColumnGroups len %d != nvec %d" % \ 3266 (len(self.groups), self.nvec)) 3267 if self.verb > verb_level: 3268 print("-- label %s: groups %s" % (label,self.groups)) 3269 elif label == 'CommandLine': 3270 self.command = data 3271 if self.verb > verb_level: 3272 print("-- label %s: command %s" % (label,self.command)) 3273 elif label == 'RowTR': 3274 self.tr = float(data) 3275 if self.verb > verb_level: 3276 print("-- label %s: TR %s" % (label,self.tr)) 3277 elif label == 'GoodList': 3278 self.goodlist = UTIL.decode_1D_ints(data) 3279 if self.goodlist: 3280 if len(self.goodlist) != self.nt: 3281 print("** GoodList missing %d rows" % \ 3282 self.nt-len(self.goodlist)) 3283 if self.verb > verb_level: 3284 print("-- label %s: goodlist %s" % (label,self.goodlist)) 3285 elif label == 'NRowFull': 3286 self.nrowfull = int(data) 3287 if self.verb > verb_level: 3288 print("-- label %s: nrowfull %s" % (label,self.nrowfull)) 3289 elif label == 'RunStart': 3290 self.runstart = UTIL.decode_1D_ints(data) 3291 if not UTIL.vals_are_sorted(self.runstart): 3292 print("** RunStart values not sorted? data = %s" % data) 3293 self.runstart = [] 3294 if self.verb > verb_level: 3295 print("-- label %s: runstart %s" % (label,self.runstart)) 3296 elif label.startswith('Basis'): 3297 self.add_basisinfo(label, data) 3298 elif self.verb > 2: 3299 print("** unknown comment label '%s'" % label) 3300 except: 3301 print("** failed to process comment label '%s'" % label) 3302 3303 if self.verb > 2 and len(self.basisinfo) > 0: 3304 for bind, bdict in enumerate(self.basisinfo): 3305 print("++ bdict %d: %s" % (bind+1, bdict)) 3306 3307 return 0 3308 3309# end Afni1D 3310# =========================================================================== 3311 3312def c1D_line2labelNdata(cline, verb=1): 3313 """expect cline to be of the form: '# LABEL = "DATA"' 3314 returning LABEL, DATA""" 3315 sline = cline.split() 3316 3317 # check for problems 3318 if len(sline) < 4 or sline[0] != "#" or sline[2] != "=": 3319 if verb > 2: print("-- skipping useless comment line: '%s'" % cline) 3320 return None, None 3321 3322 label = sline[1] # store the line label 3323 3324 dline = ' '.join(sline[3:]) # put the line back together 3325 3326 # again, check for problems 3327 if len(dline) < 2 or dline[0] != '"' or dline[-1] != '"': 3328 if verb > 1: print("** missing quotes in comment line: '%s'" % cline) 3329 return None, None 3330 3331 data = dline[1:-1] 3332 3333 if verb > 2: print("++ line2labelNdata returns '%s', '%s'" % (label,data)) 3334 3335 return label, data 3336 3337def list2_is_in_list1(list1, list2, label=''): 3338 """return 0 or 1, based on whether every element in list2 exists 3339 in list1""" 3340 3341 if not list1: return 0 3342 if not list2: return 1 3343 3344 for item in list2: 3345 if not item in list1: 3346 if label: print("-- list '%s' not subset in Afni1D" % label) 3347 return 0 3348 3349 return 1 3350 3351# end Afni1D helper functions 3352# =========================================================================== 3353 3354 3355# =========================================================================== 3356# begin AfniData - generic numeric sparse 2D float class 3357 3358g_AfniData_hist = """ 3359 30 Jul, 2010 : added AfniData class 3360 17 Aug, 2010 : get data via lib_textdata.read_data_file() 3361 - this includes modulation and duration data 3362""" 3363 3364# error constants for file tests 3365ERR_ANY_MISC = 1 # apply to errors that are not accumulated 3366ERR_ST_NEGATIVES = 2 # some times are negatives 3367ERR_ST_NON_UNIQUE = 4 # times are not unique 3368ERR_ST_NUM_RUNS = 8 # number of runs mismatch 3369ERR_ST_TOO_BIG = 16 # val >= run length 3370 3371 3372class AfniData(object): 3373 def __init__(self, filename="", mdata=None, fsl_flist=[], verb=1): 3374 """akin to a 2D float class, but do not require a square matrix 3375 3376 init from filename 3377 filename : 1D/timing file to read 3378 verb : verbose level of operations 3379 """ 3380 3381 # main variables 3382 self.fname = filename # name of data file 3383 self.name = "NoName" # more personal and touchy-feely... 3384 self.data = None # actual data (array of arrays [[]]) 3385 self.mdata = None # married data (elements are [time [mods] dur]) 3386 self.clines = None # comment lines from file 3387 self.alist = None # list of '*' counts per row 3388 3389 # descriptive variables, set from data 3390 self.nrows = 0 3391 self.ncols = 0 # non-zero if rectangular 3392 self.row_lens = [] # row lengths, if not rectangular 3393 3394 self.binary = 0 # 0/1 file? 3395 self.empty = 0 # no data at all 3396 self.married = 0 # data has modulators or durations 3397 self.mtype = 0 # married type (bits, amp, dur) 3398 self.minlen = 0 # min and max row lengths 3399 self.maxlen = 0 3400 self.dur_len = 0.0 3401 3402 self.ready = 0 # data is ready 3403 3404 # passed in variables 3405 self.tr = 0 # non-zero, if it applies 3406 self.nruns = 0 # non-zero, if known 3407 self.run_lens = [] # run lengths, in seconds or TRs 3408 self.verb = verb 3409 # self.hist = g_AfniData_hist (why did I do this?) 3410 self.write_dm = 1 # if found include durations when printing 3411 3412 # computed variables 3413 self.cormat = None # correlation mat (normed xtx) 3414 self.cosmat = None # cosine mat (scaled xtx) 3415 self.cormat_ready = 0 # correlation mat is set 3416 3417 # initialize... 3418 if self.fname: 3419 if self.init_from_filename(self.fname): return None 3420 elif mdata: 3421 if self.init_from_mdata(mdata): return None 3422 elif len(fsl_flist) > 0: 3423 if self.init_from_fsl_flist(fsl_flist): return None 3424 3425 # some accessor functions to match Afni1D 3426 def set_nruns(nruns): 3427 self.nruns = nruns 3428 3429 def set_run_lengths(run_lengths): 3430 self.row_lens = run_lengths 3431 3432 def is_rect(self): 3433 return TD.data_is_rect(self.mdata) 3434 3435 def is_square(self): 3436 if not self.is_rect(): return 0 3437 3438 rlen = len(self.mat) 3439 3440 if rlen == 0: return 1 3441 if rlen == len(self.mat[0]): return 1 3442 3443 return 0 3444 3445 def get_duration(self): 3446 """return 0 : if empty 3447 dur : if constant 3448 -1 : otherwise (existing, but not constant) 3449 """ 3450 3451 if self.mtype == 0 or len(self.mdata) == 0: return 0 3452 3453 # have actual married data 3454 3455 dur = -2 3456 for row in self.mdata: 3457 for entry in row: 3458 if dur == -2: dur = entry[2] 3459 elif dur != entry[2]: 3460 dur = -1 3461 break 3462 if dur == -1: break 3463 if dur == -2: return 0 # uninitialized means 0 3464 return dur # -1 or good value 3465 3466 def get_min_max_duration(self): 3467 """return -1, -1 : if empty 3468 min, max: otherwise 3469 """ 3470 3471 if self.mtype == 0 or len(self.mdata) == 0: return 0 3472 3473 # have actual married data 3474 3475 dmin = -1 3476 dmax = -1 3477 for row in self.mdata: 3478 for entry in row: 3479 dur = entry[2] 3480 if dmax == -1: dmin = dur # init, but not to a small thing 3481 if dur < dmin: dmin = dur 3482 if dur > dmax: dmax = dur 3483 return dmin, dmax 3484 3485 def extend_data_rows(self, newdata, promote_mtype=1, promote_rows=0): 3486 """extend each row by the corresponding row of newdata 3487 if promote_mtypes and they differ, promote so they are the same 3488 if promote_rows and they differ, promote so they are the same 3489 """ 3490 3491 if not self.ready or not newdata.ready: 3492 print('** timing elements not ready for extending rows (%d,%d)' % \ 3493 (self.ready, newdata.ready)) 3494 return 1 3495 3496 if self.nrows != newdata.nrows: 3497 print('** timing nrows differ for extending (%d, %d)' % \ 3498 (self.nrows,newdata.nrows)) 3499 if not promote_rows: return 1 3500 print(' promoting nrows between %s and %s'%(self.name, newdata.name)) 3501 self.promote_nrows(newdata) 3502 3503 if self.mtype != newdata.mtype: 3504 print('** timing elements differ in married type (%s, %s)' \ 3505 % (self.married_type_string(self.mtype), 3506 self.married_type_string(newdata.mtype))) 3507 if not promote_mtype: return 1 3508 print(' promoting mtypes between %s and %s'%(self.name, newdata.name)) 3509 self.promote_mtypes(newdata) 3510 3511 if self.verb > 1: print('++ MTiming: extending %d rows' % self.nrows) 3512 3513 for ind in range(self.nrows): 3514 if self.verb > 3: 3515 print('++ edr m #%d: extending %d rows from %d cols by %d' \ 3516 % (ind, self.nrows, len(self.data[ind]), len(newdata.data[ind]))) 3517 self.data[ind].extend(newdata.data[ind]) 3518 3519 for ind in range(self.nrows): 3520 if self.verb > 3: 3521 print('++ EDR M #%d: extending %d rows from %d cols by %d' \ 3522 % (ind, self.nrows, len(self.mdata[ind]),len(newdata.mdata[ind]))) 3523 self.mdata[ind].extend(newdata.mdata[ind]) 3524 3525 # and update ncols 3526 self.ncols = min([len(r) for r in self.data]) 3527 3528 if self.verb > 3: 3529 print(self.make_data_string(nplaces=1, flag_empty=0, check_simple=0, 3530 mesg='\nnew mdata for %s'%self.name)) 3531 print('-- min, max dur = %s %s' % self.get_min_max_duration()) 3532 3533 return 0 3534 3535 def apply_end_times_as_durs(self, newdata): 3536 """treat newdata as ending times for the events in self 3537 - sort both to be sure 3538 - whine on negatives 3539 - result will include MTYPE_DUR 3540 """ 3541 3542 if not self.ready or not newdata.ready: 3543 print('** timing elements not ready for end_times (%d,%d)' % \ 3544 (self.ready, newdata.ready)) 3545 return 1 3546 3547 if self.nrows != newdata.nrows: 3548 print('** timing nrows differ for end_times (%d, %d)' % \ 3549 (self.nrows,newdata.nrows)) 3550 3551 # check that each row is the same length 3552 okay = 1 3553 for rind in range(self.nrows): 3554 l0 = len(self.data[rind]) 3555 l1 = len(newdata.data[rind]) 3556 if l0 != l1: 3557 print('** num events differs for run %d' % (rind+1)) 3558 okay = 0 3559 if not okay: 3560 return 1 3561 3562 if self.verb > 1: print('++ MTiming: applying end_times as durs') 3563 3564 # sort both to be sure 3565 self.sort() 3566 newdata.sort() 3567 3568 # apply newdata as end times, whine if negative 3569 for rind in range(self.nrows): 3570 for eind in range(len(self.data[rind])): 3571 dur = newdata.mdata[rind][eind][0] - self.mdata[rind][eind][0] 3572 if dur < 0 and self.verb: 3573 print('** end_times as durs: have negative duration %f' % dur) 3574 self.mdata[rind][eind][2] = dur 3575 3576 # set mtype to include MTYPE_DUR (might already) 3577 self.mtype |= MTYPE_DUR 3578 3579 return 0 3580 3581 def promote_nrows(self, newdata, dataval=None, maryval=None): 3582 """if differing nrows, extend to match 3583 3584 if dataval == None: new rows are empty 3585 else: pad by ncols vector of dataval elements 3586 """ 3587 3588 diff = self.nrows - newdata.nrows 3589 3590 if diff == 0: return 0 3591 3592 if dataval == None: df = [] 3593 else: df = [dataval for i in range(self.ncols)] 3594 if maryval == None: mf = [] 3595 else: mf = [maryval for i in range(self.ncols)] 3596 3597 if diff < 0: updater = self 3598 else : updater = newdata 3599 3600 for r in range(abs(diff)): 3601 updater.data.append(df) 3602 updater.mdata.append(mf) 3603 3604 return 0 3605 3606 def promote_mtypes(self, newdata): 3607 """if married types differ, promote data to match""" 3608 3609 # if they are the same, there is nothing to do 3610 if self.mtype == newdata.mtype: return 0 3611 3612 # else they will certainly be married after this, 3613 # so set both to married, and mtype to the combined one 3614 self.married = 1 3615 newdata.married = 1 3616 new_mtype = self.mtype | newdata.mtype 3617 self.mtype = new_mtype 3618 newdata.mtype = new_mtype 3619 3620 if self.verb > 2: 3621 print('++ promoting married type to %d (%s from %s)' \ 3622 % (new_mtype, self.name, newdata.name)) 3623 3624 # try to find some data 3625 sval = [] 3626 for v in self.mdata: 3627 if len(v) > 0: 3628 sval = v[0] 3629 break 3630 nval = [] 3631 for v in newdata.mdata: 3632 if len(v) > 0: 3633 nval = v[0] 3634 break 3635 3636 # if either has no data, we are done 3637 if len(sval) == 0 or len(nval) == 0: return 0 3638 3639 # so we have data, see who has more amplitude modulators 3640 3641 adiff = len(sval[1]) - len(nval[1]) 3642 if adiff < 0: self.pad_amp_mods(-adiff) 3643 elif adiff > 0: newdata.pad_amp_mods(adiff) 3644 if adiff != 0: # promote mtypes 3645 self.mtype = (self.mtype | MTYPE_AMP) 3646 newdata.mtype = (newdata.mtype | MTYPE_AMP) 3647 3648 # nothing to do for MTYPE_DUR, except promote mtype, done above 3649 3650 return 0 3651 3652 def pad_amp_mods(self, ext_len): 3653 """extend the amplitudes by length ext_len list of 1""" 3654 if ext_len <= 0: return 3655 if self.verb > 2: 3656 print('++ pad_amp_mods(%d) for %s' % (ext_len, self.name)) 3657 elist = [1] * ext_len 3658 for row in self.mdata: 3659 for val in row: 3660 val[1].extend(elist) 3661 3662 def show_married_info(self): 3663 print('-- modulation type: %s' % self.married_info_string()) 3664 3665 def married_info_string(self): 3666 if self.mtype == MTYPE_NONE: return 'not a married format' 3667 3668 namp = self.num_amplitudes() 3669 adur = self.ave_dur_modulation() 3670 3671 if self.mtype == MTYPE_AMP: 3672 return 'amplitude modulation (%d modulators)' % namp 3673 3674 if self.mtype == MTYPE_DUR: 3675 return 'duration modulation (average duration %g)' % adur 3676 3677 if self.mtype == MTYPE_AMP|MTYPE_DUR: 3678 return 'amp and dur modulation (%d amp mods, ave dur %g)'%(namp, adur) 3679 3680 return '** invalid modulation type %d' % self.mtype 3681 3682 def num_amplitudes(self): 3683 if not self.mtype & MTYPE_AMP: return 0 3684 if not self.mdata: return 0 3685 for row in self.mdata: 3686 if len(row) == 0: continue 3687 return len(row[0][1]) # have amplitudes, return the length 3688 return 0 # no valid rows found 3689 3690 def ave_dur_modulation(self): 3691 if not self.mtype & MTYPE_DUR: return 0 3692 if not self.mdata: return 0 3693 lsum, count = 0.0, 0 3694 for row in self.mdata: 3695 if len(row) == 0: continue 3696 count += len(row) 3697 lsum += UTIL.loc_sum([entry[2] for entry in row]) 3698 if count == 0: return 0 3699 return lsum*1.0/count # be wary of integers 3700 3701 def married_type_string(self, mtype=None): 3702 if mtype == None: mtype = self.mtype 3703 3704 if mtype == MTYPE_NONE: return 'None' 3705 elif mtype == MTYPE_AMP: return 'Amp Mod' 3706 elif mtype == MTYPE_DUR: return 'Dur Mod' 3707 elif mtype == MTYPE_AMP | MTYPE_DUR: return 'Amp/Dur Mod' 3708 else: return 'Unknown' 3709 3710 def randomize_trs(self, seed=0): 3711 """reorder the matrix rows randomly""" 3712 if self.verb > 1: print('-- randomizing %d trs (seed=%d)' \ 3713 % (len(self.data, seed))) 3714 if seed: 3715 import random 3716 random.seed(seed) 3717 UTIL.shuffle(self.data) 3718 3719 def sort(self, rev=0): 3720 """sort each row (optionally reverse order)""" 3721 if not self.ready: return 1 3722 3723 if self.verb > 1: print('-- sorting AfniData ...') 3724 3725 for row in self.data: 3726 if rev: row.sort(reverse=True) 3727 else: row.sort() 3728 3729 for row in self.mdata: 3730 if rev: row.sort(reverse=True) 3731 else: row.sort() 3732 3733 return 0 3734 3735 def transpose(self): 3736 """the tranpose operation requires rectangular data""" 3737 if not self.ready: return 1 3738 3739 # allow transpose if max row length is 1 3740 if self.maxlen > 1 and not self.is_rect(): 3741 print('** cannot take transpose, data is not rectangular') 3742 return 1 3743 3744 if self.verb > 1: print('-- AData: taking transpose...') 3745 3746 if self.empty: return 0 3747 3748 # get mdata and data 3749 3750 # if column (maybe non-rectagular), allow transpose 3751 # (akin to when 3dDeconvolve might mis-interpret local times as global) 3752 if self.maxlen == 1 and not self.is_rect(): 3753 newdata = [] 3754 for row in self.data: newdata.extend(row) 3755 self.data = [newdata] 3756 newdata = [] 3757 for row in self.mdata: newdata.extend(row) 3758 self.mdata = [newdata] 3759 else: # typical rectagular case 3760 newdata = [] 3761 for col in range(len(self.data[0])): 3762 newdata.append([self.data[row][col] for row in range(self.nrows)]) 3763 del(self.data) 3764 self.data = newdata 3765 3766 newdata = [] 3767 for col in range(len(self.mdata[0])): 3768 newdata.append([self.mdata[row][col] for row in range(self.nrows)]) 3769 del(self.mdata) 3770 self.mdata = newdata 3771 3772 self.nrows = len(self.mdata) 3773 self.ncols = len(self.mdata[0]) 3774 3775 return 0 3776 3777 def copy(self): 3778 """return a complete (deep)copy of the current AfniTiming""" 3779 return copy.deepcopy(self) 3780 3781 def is_empty(self): 3782 """empty is either nrows == 0 or each row is empty""" 3783 if not self.ready: return 1 # seems safer 3784 3785 if self.nrows < 1: return 1 3786 3787 # now check each row 3788 for row in self.mdata: 3789 if len(row) > 0: return 0 # found something 3790 3791 return 1 3792 3793 def make_single_row_string(self, row=-1, nplaces=3, mplaces=-1, 3794 flag_empty=0, simple=1): 3795 """return a string of row data, to the given number of decimal places 3796 if row is non-negative, return a string for the given row 3797 3798 if not simple, show_durs controls display of durations 3799 """ 3800 if not self.ready: return '' 3801 if row < 0 or row >= self.nrows: 3802 if self.verb > 0: print('** row %d out of range for printing' % row) 3803 return '' 3804 3805 data = self.mdata[row] 3806 rstr = '' 3807 if self.verb > 2 and not flag_empty: rstr += 'run %02d : ' % (row+1) 3808 3809 # if flagging an empty run, use '*' characters 3810 # (little gain in trying to usually put one '*') 3811 if len(data) == 0 and flag_empty: rstr += '* *' 3812 3813 for val in data: 3814 if simple: 3815 if nplaces >= 0: rstr += '%.*f ' % (nplaces, val[0]) 3816 else: rstr += '%g ' % (val[0]) 3817 else: 3818 if self.mtype & MTYPE_AMP and len(val[1]) > 0: 3819 if mplaces >= 0: alist = ['*%.*f'%(nplaces, a) for a in val[1]] 3820 else: alist = ['*%g'%a for a in val[1]] 3821 astr = ''.join(alist) 3822 else: astr = '' 3823 3824 # if married and want durations, include them 3825 if self.write_dm and (self.mtype & MTYPE_DUR): 3826 if mplaces >= 0: dstr = ':%.*f' % (nplaces, val[2]) 3827 else: dstr = ':%g' % val[2] 3828 else: dstr = '' 3829 3830 if nplaces >= 0: rstr += '%.*f%s%s ' % (nplaces, val[0], astr, dstr) 3831 else: rstr += '%g%s%s ' % (val[0], astr, dstr) 3832 3833 return rstr + '\n' 3834 3835 def make_data_string(self, row=-1, nplaces=3, mplaces=-1, 3836 flag_empty=0, check_simple=1, mesg=''): 3837 """return a string of row data, to the given number of decimal places 3838 if row is non-negative, return a string for the given row, else 3839 return a string of all rows 3840 row : make a string for just the single row 3841 nplaces : number of decimal places to show 3842 mplaces : number of decimal places to show in married fields 3843 flag_empty : if empty row, use the '*' format 3844 check_simple : if set and dur_len=0, use timing format 3845 mesg : display the message before data 3846 """ 3847 3848 if check_simple and self.dur_len == 0.0: simple = 1 3849 elif self.mtype == 0: simple = 1 3850 else: simple = 0 3851 3852 # init return string based on message 3853 if len(mesg) > 0: rstr = "%s :\n" % mesg 3854 else: rstr = '' 3855 3856 if self.verb > 2: 3857 print('-- writing %s, simple=%d, wdm=%d, const=%d, dur_len=%d' \ 3858 % (self.name, simple, self.write_dm, self.durs_are_constant(), 3859 self.dur_len)) 3860 3861 if row >=0: 3862 return rstr+self.make_single_row_string(row, nplaces=nplaces, 3863 mplaces=mplaces, flag_empty=flag_empty, simple=simple) 3864 3865 # make it for all rows 3866 for ind in range(self.nrows): 3867 rstr += self.make_single_row_string(ind, nplaces=nplaces, 3868 mplaces=mplaces, flag_empty=flag_empty, simple=simple) 3869 3870 return rstr 3871 3872 def durs_are_constant(self): 3873 """just check mdata""" 3874 if self.mdata == None: return 1 3875 dall = [] 3876 for mrow in self.mdata: 3877 dall.extend([e[2] for e in mrow]) 3878 return UTIL.vals_are_constant(dall) 3879 3880 def check_constant_duration(self): 3881 """return whether constant, and duration if so 3882 (akin to durs_are_constant(), but also return that duration) 3883 return -1 for duration, if unknown 3884 """ 3885 cdur = -1 3886 if self.mdata == None: return 1, cdur 3887 dall = [] 3888 for mrow in self.mdata: 3889 dall.extend([e[2] for e in mrow]) 3890 3891 if UTIL.vals_are_constant(dall): 3892 if len(dall) > 0: cdur = dall[0] 3893 rv = 1 3894 else: 3895 rv = 0 3896 3897 return rv, cdur 3898 3899 def write_as_timing(self, fname='', nplaces=-1, mplaces=-1, check_simple=1): 3900 """write the current M timing out, with nplaces right of the decimal""" 3901 if not self.ready: 3902 print('** M Timing: not ready to write') 3903 return 1 3904 3905 if fname == '': fname = self.fname 3906 3907 if fname == '-' or fname == 'stdout': 3908 fp = sys.stdout 3909 else: 3910 fp = open(fname, 'w') 3911 if not fp: 3912 print("** failed to open '%s' for writing Mtiming" % fname) 3913 return 1 3914 3915 if self.verb > 1: 3916 print("++ writing %d MTiming rows to %s" % (self.nrows, fname)) 3917 3918 fp.write(self.make_data_string(nplaces=nplaces, mplaces=mplaces, 3919 flag_empty=1, check_simple=check_simple)) 3920 3921 if fp != sys.stdout: 3922 fp.close() 3923 3924 return 0 3925 3926 def looks_like_1D(self, run_lens=[], nstim=0, verb=1): 3927 """check if the file looks like 1D 3928 - if so and verb, show warnings, too""" 3929 3930 # negative result is terminal, positive can continue with warnings 3931 errs = self.file_type_errors_1D(run_lens=run_lens, nstim=nstim, verb=verb) 3932 if errs < 0: 3933 if verb>0: print('== BAD: %s does not look like valid 1D'%self.fname) 3934 return 0 3935 else: 3936 if verb > 0: 3937 self.file_type_warnings_1D(run_lens=run_lens, nstim=nstim) 3938 print('== GOOD: %s looks like 1D' % self.fname) 3939 if errs > 0: return 0 3940 return 1 3941 3942 3943 def file_type_errors_1D(self, run_lens=[], nstim=0, verb=1): 3944 """return whether errors exist for 1D (stim_file) format 3945 - data should be rectangular 3946 (what else can we even check?) 3947 if run_lens is passed, 3948 - nrows should be at least sum(run_lens) 3949 if nstim is passed, 3950 - number of columns should equal nstim 3951 3952 return -1 on fatal error 3953 0 on OK 3954 1 on non-fatal error 3955 """ 3956 3957 if not self.ready: 3958 print('** looks_like_1D: data not ready') 3959 return -1 3960 3961 if self.empty: 3962 if verb > 1: print("-- empty file %s okay as 1D file" % self.fname) 3963 return 0 3964 3965 # error occur only as singles, to test against verb > 1 3966 3967 errors = 0 3968 if not self.is_rect(): 3969 errors |= ERR_ANY_MISC 3970 if verb > 1: print("** file %s is not rectangular" % self.fname) 3971 3972 # keep same local variable 3973 rlens = run_lens 3974 3975 nruns = len(rlens) 3976 if nruns > 0: 3977 # if TR, scale it in 3978 tot_dur = UTIL.loc_sum(rlens) 3979 3980 # if nrows does not match total run length, fail 3981 if tot_dur != self.nrows: 3982 errors |= ERR_ANY_MISC 3983 if verb > 1: 3984 print("** file %s: nrows does not match run dur: %d != %d" \ 3985 % (self.fname, self.nrows, tot_dur)) 3986 3987 if nstim > 0 and nstim != self.ncols: 3988 errors |= ERR_ANY_MISC 3989 if verb > 1: print("** file %s: ncols %d != nstim %d" \ 3990 % (self.fname, self.ncols, nstim)) 3991 3992 if errors: return -1 3993 else: return 0 3994 3995 def file_type_warnings_1D(self, run_lens=[], nstim=0): 3996 """akin to ft_errors_1D, but just show any warnings 3997 if run_lens is passed, 3998 - warn if too many 3999 if nstim is passed, 4000 - number of columns should equal nstim 4001 return whether warnings are shown 4002 """ 4003 4004 if self.file_type_errors_1D(run_lens=run_lens, nstim=nstim, verb=0) < 0: 4005 print('** file %s is not valid as a 1D format' % self.name) 4006 return 1 4007 4008 if self.empty: 4009 print("-- empty file %s okay as 1D file" % self.fname) 4010 return 0 4011 4012 # keep same local variable 4013 rlens = run_lens 4014 4015 nruns = len(rlens) 4016 if nruns > 0: 4017 # if TR, scale it in 4018 tot_dur = UTIL.loc_sum(rlens) 4019 4020 if tot_dur != self.nrows: 4021 print("** warning for 1D file %s, num rows != num TRs: %d != %d" \ 4022 % (self.fname, self.nrows, tot_dur)) 4023 return 1 4024 4025 return 0 4026 4027 def looks_like_local_times(self, run_lens=[], tr=0.0, verb=1): 4028 """check if the file looks like local times 4029 - if so and verb, show warnings, too""" 4030 4031 # negative result is terminal, positive can continue with warnings 4032 errs = self.file_type_errors_local(run_lens=run_lens, tr=tr, verb=verb) 4033 if errs < 0: 4034 if verb > 0: 4035 print('== BAD: %s does not look like local stim_times' % self.fname) 4036 return 0 4037 else: 4038 if verb > 0: 4039 self.file_type_warnings_local(run_lens=run_lens, tr=tr) 4040 print('== GOOD: %s looks like local stim_times' % self.fname) 4041 if errs > 0: return 0 4042 return 1 4043 4044 def file_type_errors_local(self, run_lens=[], tr=0.0, verb=1): 4045 """return whether data has errors for local stim_times format 4046 - times should be non-negative 4047 - times should be unique per run 4048 if run_lens is passed, 4049 - number of runs should match nrows 4050 - maximum should be less than current run_length 4051 if tr is passed, scale the run lengths 4052 4053 return -1 on fatal erros 4054 0 on OK 4055 1 on non-fatal errors 4056 """ 4057 4058 # possibly scale run_lengths 4059 if tr > 0.0: rlens = [tr*rl for rl in run_lens] 4060 else: rlens = run_lens 4061 4062 nruns = len(rlens) 4063 4064 if not self.ready: 4065 print('** looks_like_local_times: data not ready') 4066 return -1 4067 4068 if self.empty: 4069 if verb > 1: print("-- empty file %s okay as local_times" % self.fname) 4070 return 0 4071 4072 # in case we know nothing, just check that values are non-negative 4073 # and unique 4074 # - if we have run lengths, check the number of runs and maximums, too 4075 errors = 0 4076 4077 # allow bad nruns (but warn) 4078 4079 for rind in range(len(self.data)): 4080 # start with row copy 4081 row = self.data[rind][:] 4082 if len(row) == 0: continue 4083 row.sort() 4084 first = row[0] 4085 last = row[-1] 4086 # allow negatives (but warn) 4087 if not UTIL.vals_are_increasing(row): 4088 errors |= ERR_ST_NON_UNIQUE 4089 if verb > 1: print("** file %s, row %d has repetitive times" \ 4090 % (self.fname, rind)) 4091 # allow times after end of run (but warn) 4092 4093 del(row) 4094 4095 # no fatal errors yet 4096 if errors: return 1 4097 else: return 0 4098 4099 def looks_local_but_3dD_global(self, warn=0, maxruns=20): 4100 """return 1 if the timing looks local (has '*' anywhere), 4101 but would be read as global by 3dDeconvolve (max cols == 1) 4102 4103 true if '*' entries exist, and: 4104 - max entries (per row) is 1 (including '*') 4105 - there is some time after row 1 4106 4107 true if '*' entries do not exist, and: 4108 - there is exactly one time per row 4109 - either of: 4110 - len(run_lens) == nrows 4111 - there are <= maxruns rows 4112 4113 if verb: warn user 4114 """ 4115 4116 if not self.ready: 4117 print('** LLB3G: data not ready') 4118 return 0 4119 4120 if self.empty: return 0 4121 4122 # has_alist will imply it exists, is valid and there are '*' entries 4123 has_alist = (self.alist != None) 4124 if has_alist: has_alist = (len(self.alist) == len(self.data)) 4125 if has_alist: has_alist = (UTIL.loc_sum(self.alist) > 0) 4126 4127 minntimes = 2 # at which point we would not care 4128 maxent = 0 4129 late_times = 0 4130 4131 # track max entries per row, including '*' chars 4132 for rind, row in enumerate(self.data): 4133 nent = len(row) 4134 if rind > 0 and nent > 0: late_times = 1 # stim exist after row 0 4135 if nent < minntimes: minntimes = nent # min stim times per row 4136 if has_alist: nent += self.alist[rind] # total events per row 4137 if nent > maxent: maxent = nent # max events across rows 4138 4139 # save a little time if we are safe 4140 if maxent > 1: return 0 4141 4142 # if maxent is not 1, there is no problem 4143 if maxent != 1: return 0 4144 4145 # if there is no event after row 1, there is no problem 4146 if not late_times: return 0 4147 4148 # we have at most 1 entry per row, and times after row 1, 4149 # get to the stated tests 4150 4151 if has_alist: 4152 if warn: 4153 print("** timing file %s looks like local times from '*', but\n" \ 4154 " might be interpreted as global times by 3dDeconvovle\n"\ 4155 " because it has only one column\n" \ 4156 " (consider adding one '*', giving that row 2 entries)\n"\ 4157 % self.fname) 4158 return 1 4159 4160 if len(self.run_lens) > 0: 4161 if len(self.run_lens) == self.nrows: 4162 if warn: 4163 print("** timing file %s looks like global times, but\n" \ 4164 " Nruns == Nstim, so maybe it is local\n" \ 4165 " (if local, add '*' to some row)\n" \ 4166 % self.fname) 4167 return 1 4168 elif self.nrows > 0 and self.nrows < maxruns: 4169 if warn: 4170 print("** timing file %s looks like global times, but\n" \ 4171 " has very few stim, so might be local\n" \ 4172 " (if local, add '*' to some row)\n" \ 4173 % self.fname) 4174 return 1 4175 4176 return 0 4177 4178 def file_type_warnings_local(self, run_lens=[], tr=0.0): 4179 """warn about any oddities in local timing 4180 - times should be non-negative 4181 - times should be unique per run 4182 if run_lens is passed, 4183 - number of runs should match nrows 4184 - maximum should be less than current run_length 4185 if tr is passed, scale the run lengths 4186 return whether any warnings are printed 4187 """ 4188 4189 # possibly scale run_lengths 4190 if tr > 0.0: rlens = [tr*rl for rl in run_lens] 4191 else: rlens = run_lens 4192 4193 nruns = len(rlens) 4194 4195 if not self.ready: 4196 print('** looks_like_local_times_warn: data not ready') 4197 return 1 4198 4199 if self.file_type_errors_local(run_lens=run_lens,tr=tr,verb=0) < 0: 4200 print('** file %s is not valid as local times format' % self.name) 4201 return 1 4202 4203 if self.empty: 4204 print("-- empty file %s okay as local_times" % self.fname) 4205 return 0 4206 4207 # make a list of warnings to print at the end 4208 warnings = [] 4209 4210 # - if we have run lengths, check the number of runs and maximums, too 4211 if nruns > 0: 4212 if nruns != self.nrows: 4213 warnings.append(" - %d rows does not match %d runs" \ 4214 % (self.nrows, nruns)) 4215 4216 check_alist = (self.alist != None) 4217 if check_alist: check_alist = len(self.alist) == len(self.data) 4218 maxent = 0 4219 4220 wcount = [0,0,0] # limit warnings of each type (to 2 for now) 4221 for rind in range(len(self.data)): 4222 # start with row copy 4223 row = self.data[rind][:] 4224 4225 # track max entries per row, including '*' chars 4226 rlen = len(row) 4227 nent = rlen 4228 if check_alist: nent += self.alist[rind] 4229 if nent > maxent: maxent = nent 4230 4231 if rlen == 0: continue 4232 4233 row.sort() 4234 first = row[0] 4235 last = row[-1] 4236 if first < 0: 4237 wcount[0] += 1 4238 if wcount[0] <= 2: 4239 warnings.append(" - row %d has negative time %g"%(rind, first)) 4240 if not UTIL.vals_are_increasing(row): 4241 wcount[1] += 1 4242 if wcount[1] <= 2: 4243 warnings.append(" - row %d has repetitive times" % (rind)) 4244 if nruns > 0 and rind < nruns: 4245 if last >= rlens[rind]: 4246 wcount[2] += 1 4247 if wcount[2] <= 2: 4248 warnings.append(" - row %d : time %g exceeds run dur %g" \ 4249 % (rind, last, rlens[rind])) 4250 del(row) 4251 4252 # if any type exceeded the limit (of 2), print a general count 4253 if wcount[0] > 2: 4254 warnings.append(" * %d rows with negative times ..."%wcount[0]) 4255 if wcount[1] > 2: 4256 warnings.append(" * %d rows with repetitive times ..."%wcount[1]) 4257 if wcount[2] > 2: 4258 warnings.append(" * %d row times exceed run dur %g ..." \ 4259 % (wcount[2], rlens[rind])) 4260 4261 if maxent == 1 and self.nrows > 1: 4262 awarn = 0 4263 if check_alist: 4264 if UTIL.loc_sum(self.alist) > 0: awarn = 1 4265 if awarn: warnings.append( \ 4266 " * single column looks local from '*',\n" \ 4267 " but 3dDeconvolve would interpret as global") 4268 else: warnings.append( \ 4269 " * 3dDeconvolve would interpret single column as global") 4270 4271 if len(warnings) > 0: 4272 print('** warnings for local stim_times format of file %s' % self.fname) 4273 for w in warnings: print(w) 4274 return 1 4275 else: return 0 4276 4277 4278 def looks_like_global_times(self, run_lens=[], tr=0.0, verb=1): 4279 """check if the file looks like global times 4280 - if so and verb, show warnings, too""" 4281 4282 # negative result is terminal, positive can continue with warnings 4283 errs = self.file_type_errors_global(run_lens=run_lens, tr=tr, verb=verb) 4284 if errs < 0: 4285 if verb > 0: 4286 print('== BAD: %s does not look like global stim_times' % self.fname) 4287 return 0 4288 else: 4289 if verb > 0: 4290 self.file_type_warnings_global(run_lens=run_lens, tr=tr) 4291 print('== GOOD: %s looks like global stim_times' % self.fname) 4292 if errs > 0: return 0 4293 return 1 4294 4295 def file_type_errors_global(self, run_lens=[], tr=0.0, verb=1): 4296 """ return -1 on fatal erros 4297 0 on OK 4298 1 on non-fatal errors 4299 """ 4300 4301 if not self.ready: 4302 print('** looks_like_1D: data not ready') 4303 return -1 4304 4305 if self.empty: 4306 if verb > 1: print("-- empty file %s okay as global_times" % self.fname) 4307 return 0 4308 4309 errors = 0 4310 ferrors = 0 4311 4312 # must be one row or column 4313 if self.ncols != 1 and self.nrows != 1: 4314 errors |= ERR_ANY_MISC 4315 if verb > 1: print("** file %s is not a single column" % self.fname) 4316 4317 # must be rectangular 4318 if not self.is_rect(): 4319 errors |= ERR_ANY_MISC 4320 if verb > 1: print("** file %s is not a rectangular" % self.fname) 4321 4322 # this is fatal, as opposed to nrows/ncols 4323 if self.maxlen > 1: 4324 ferrors |= ERR_ANY_MISC 4325 if verb: print('** file %s has rows longer than 1' % self.fname) 4326 4327 4328 # negative times are not errors, but warnings 4329 4330 # possibly scale run_lengths 4331 4332 # note the total duration (tr == -1 implies just count TRs) 4333 4334 # late times are not errors, but warnings 4335 4336 # repetition times are not errors, but warnings 4337 4338 if ferrors: return -1 4339 elif errors: return 1 4340 else: return 0 4341 4342 def file_type_warnings_global(self, run_lens=[], tr=0.0, verb=1): 4343 4344 if not self.ready: 4345 print('** looks_like_1D: data not ready') 4346 return 1 4347 4348 if self.empty: 4349 if verb > 1: print("-- empty file %s okay as global_times" % self.fname) 4350 return 0 4351 4352 if self.file_type_errors_global(run_lens=run_lens,tr=tr,verb=0) < 0: 4353 print('** file %s is not valid as global times format' % self.name) 4354 return 1 4355 4356 # make a list of warnings to print at the end 4357 warnings = [] 4358 4359 # get a single sequence of numbers, depending on the direction 4360 if self.nrows == 1: data = self.data[0] 4361 else:data = [row[0] for row in self.data if len(row)>0] 4362 4363 data.sort() 4364 4365 if data[0] < 0.0: 4366 warnings.append(" - negative stim time %g" % (data[0])) 4367 4368 # possibly scale run_lengths 4369 if tr > 0.0: rlens = [tr*rl for rl in run_lens] 4370 else: rlens = run_lens 4371 4372 # note the total duration (tr == -1 implies just count TRs) 4373 endoftime = UTIL.loc_sum(rlens) 4374 if tr > 0.0: endoftime *= tr 4375 4376 if data[-1] >= endoftime: 4377 warnings.append(" - time %g after all runs, %g" 4378 % (data[-1], endoftime)) 4379 4380 if not UTIL.vals_are_increasing(data): 4381 warnings.append(" - has repeated times") 4382 4383 if len(warnings) > 0: 4384 print('** warnings for global stim_times format of file %s'%self.fname) 4385 for w in warnings: print(w) 4386 return 1 4387 else: return 0 4388 4389 def init_from_filename(self, fname): 4390 """file could be 1D, timing or married timing data 4391 4392 For now, store complete result but focus on times only. 4393 """ 4394 4395 mdata, clines, alist = TD.read_married_file(fname, verb=self.verb) 4396 if mdata == None: 4397 if self.verb > 0: print('** A1D: failed to read data file %s' % fname) 4398 return 1 4399 4400 # init name from filename 4401 aname = BASE.afni_name(self.fname) 4402 self.name = aname.pve() 4403 4404 # note whether the data is married (modulation or duration) 4405 self.mtype = TD.married_type(mdata) 4406 if self.mtype: self.married = 1 4407 4408 # data will ignore any married information 4409 self.data = [[val[0] for val in row] for row in mdata] 4410 self.mdata = mdata 4411 self.clines = clines 4412 self.fname = fname 4413 self.alist = alist 4414 4415 self.nrows = len(self.data) 4416 self.row_lens = [len(row) for row in self.data] 4417 4418 # empty data includes existing but empty runs 4419 if len(self.data) == 0: 4420 self.maxlen = 0 4421 self.minlen = 0 4422 else: 4423 self.maxlen = max([len(drow) for drow in self.data]) 4424 self.minlen = min([len(drow) for drow in self.data]) 4425 4426 # accept an empty file? 4427 if self.nrows == 0 or self.maxlen == 0: 4428 self.empty = 1 4429 self.ready = 1 4430 return 0 4431 4432 # if row lengths are all the same, use ncols, instead 4433 if UTIL.vals_are_constant(self.row_lens): 4434 self.ncols = self.row_lens[0] 4435 del(self.row_lens) 4436 self.row_lens = [] 4437 4438 # check to see if it is a 0/1 file 4439 self.binary = 1 4440 for row in self.data: 4441 if not UTIL.vals_are_0_1(row): 4442 self.binary = 0 4443 break 4444 4445 self.ready = 1 4446 4447 return 0 4448 4449 def init_from_fsl_flist(self, flist): 4450 """Convert FSL flist files into a married data structure, 4451 and then init from that. 4452 4453 Files in flist are considered to be in the FSL format: 4454 time duration amplitude 4455 and are considered to be one run per file. 4456 4457 --> So basically just swap amp and dur, and put amp in list. 4458 4459 Output is a array of mdata lists, of the form: 4460 [ 4461 [ [time [AM list] duration] [] [] ... ] 4462 [ [time [AM list] duration] [] [] ... ] 4463 ] 4464 4465 mdata[0] is the event list for the first run 4466 mdata[0][0] is the first event for the first run 4467 mdata[0][0][0] is the first event start time 4468 mdata[0][0][1] is the array of first event amplitudes 4469 mdata[0][0][1][0] is the first amplitude of the first event 4470 mdata[0][0][2] is the first event duration 4471 4472 FSL files have only single amplitudes 4473 """ 4474 4475 if len(flist) < 1: return 1 4476 4477 mdata = [] 4478 amp_all = [] # track all amplitudes 4479 4480 if self.verb > 1: 4481 print('-- reading FSL timing for %d runs' % len(flist)) 4482 4483 for find, fname in enumerate(flist): 4484 try: td = TD.read_1D_file(fname) 4485 except: 4486 print("** failed to read FSL timing file %d, '%s'" % (find, fname)) 4487 return 1 4488 if td == None: 4489 print("** failed to read FSL timing file %d, '%s'" % (find, fname)) 4490 return 1 4491 4492 # check for an empty run, either nothing in the file, 4493 # or the special case of 0, 0, 0 4494 if len(td) == 0: 4495 if self.verb > 2: 4496 print('-- FSL timing file %s has no events' % fname) 4497 mdata.append([]) 4498 continue 4499 elif len(td) == 1: 4500 if UTIL.vals_are_constant(td[0], cval=0): 4501 if self.verb > 2: 4502 print('-- FSL timing file %s shows 1 empty event' % fname) 4503 mdata.append([]) 4504 continue 4505 4506 # 1 entry is time, 2 includes duration, 3 includes amplitude 4507 if len(td[0]) == 0: 4508 if self.verb > 2: 4509 print('-- FSL timing file %s missing events' % fname) 4510 elist = [] 4511 elif len(td[0]) == 1: 4512 elist = [[ev[0], [1], 0] for ev in td] 4513 elif len(td[0]) == 2: 4514 elist = [[ev[0], [1], ev[1]] for ev in td] 4515 else: 4516 elist = [[ev[0], [ev[2]], ev[1]] for ev in td] 4517 # and track all modulators 4518 amp_all.extend([ev[2] for ev in td]) 4519 4520 mdata.append(elist) 4521 4522 # if all amplitudes are constand 0 or 1, remove them 4523 if UTIL.vals_are_constant(amp_all, cval=0): cval = 0 4524 elif UTIL.vals_are_constant(amp_all, cval=1): cval = 1 4525 else: cval = 2 4526 if cval <= 1: 4527 if self.verb > 1: 4528 print('-- FSL amplitudes are all %s, removing...' % cval) 4529 for mrow in mdata: 4530 for event in mrow: 4531 event[1] = [] 4532 4533 return self.init_from_mdata(mdata) 4534 4535 def init_from_mdata(self, mdata): 4536 """mdata should be of the form: 4537 one row per run, where each row is a list of events: 4538 [time [AM list] duration] 4539 """ 4540 4541 if not self.mdata_looks_valid(mdata): 4542 return 1 4543 4544 self.mdata = mdata 4545 4546 # note whether the data is married (modulation or duration) 4547 self.mtype = TD.married_type(mdata) 4548 if self.mtype: self.married = 1 4549 4550 # if durs, but constant, set dur_len instead 4551 if self.mtype & MTYPE_DUR: 4552 isconst, dlen = self.check_constant_duration() 4553 if isconst: 4554 self.write_dm = 0 4555 self.dur_len = dlen 4556 else: 4557 self.dur_len = -1 4558 4559 # data will ignore any married information 4560 self.data = [[val[0] for val in row] for row in mdata] 4561 self.clines = None 4562 4563 # init alist to be 0, 1 or 2, for each run so at least 2 "events" 4564 self.alist = [0] * len(mdata) 4565 for rind, run in enumerate(mdata): 4566 if len(run) == 0: 4567 self.alist[rind] = 2 4568 elif len(run) == 1: 4569 self.alist[rind] = 1 4570 4571 self.nrows = len(self.data) 4572 self.row_lens = [len(row) for row in self.data] 4573 4574 # empty data includes existing but empty runs 4575 if len(self.data) == 0: 4576 self.maxlen = 0 4577 self.minlen = 0 4578 else: 4579 self.maxlen = max([len(drow) for drow in self.data]) 4580 self.minlen = min([len(drow) for drow in self.data]) 4581 4582 # accept an empty file? 4583 if self.nrows == 0 or self.maxlen == 0: 4584 self.empty = 1 4585 self.ready = 1 4586 return 0 4587 4588 # if row lengths are all the same, use ncols, instead 4589 if UTIL.vals_are_constant(self.row_lens): 4590 self.ncols = self.row_lens[0] 4591 del(self.row_lens) 4592 self.row_lens = [] 4593 4594 # check to see if it is a 0/1 file 4595 self.binary = 1 4596 for row in self.data: 4597 if not UTIL.vals_are_0_1(row): 4598 self.binary = 0 4599 break 4600 4601 self.ready = 1 4602 4603 return 0 4604 4605 def mdata_looks_valid(self, mdata, verb=0): 4606 """show be rows of [float [float list] float]""" 4607 errs = 0 4608 for rind, row in enumerate(mdata): 4609 for eind, event in enumerate(row): 4610 # check for bad event 4611 if len(event) != 3: 4612 if verb <= 0: return 0 # quiet failure 4613 errs += 1 4614 print('** (0-based) run %d, event %d: bad form: %s' \ 4615 % (rind, eind, event)) 4616 if verb < 2: return 0 4617 if errs: return 0 4618 return 1 4619 4620 def select_runs(self, rlist): 4621 """modify the timing element where the new one is based off the 4622 listed runs of the old one 4623 4624 Elements of rlist are 1-based run indices from the old timing, 4625 itemizing the runs of the new timing. A run index of zero means 4626 to make it empty. 4627 4628 So rlist should contain integers in [0, old_num_runs]. 4629 4630 return 0 on success 4631 """ 4632 4633 nnew = len(rlist) 4634 4635 if not self.ready: 4636 print("** select_runs: not ready") 4637 return 1 4638 4639 if nnew < 1: 4640 print("** select_runs: missing run inputs (use 0 for empty)") 4641 return 1 4642 4643 # -------------------------------------------------- 4644 # limit rlist elements to [0, nold] 4645 # (but allow an index of zero to mean an empty run) 4646 4647 # note the limit on rlist entries 4648 nold = 0 4649 if self.data != None: 4650 nold = len(self.data) 4651 4652 if max(rlist) > nold or min(rlist) < 0: 4653 print("** select_runs, bad rlist limits in: %s" % rlist) 4654 return 1 4655 4656 # update number of runs to match new list 4657 self.nruns = nnew 4658 self.nrows = nnew 4659 4660 # -------------------------------------------------- 4661 # then just perform the same operation on every list 4662 # (base variable should appear at top and bottom) 4663 # (else:append() value will depend on the list) 4664 4665 # self.data 4666 lorig = self.data 4667 if lorig != None: 4668 d = [] 4669 if self.verb > 3: print("++ padding data") 4670 for rind in rlist: 4671 # if > 0, get a full copy of an old run, else use an empty one 4672 if rind > 0: d.append(lorig[rind-1][:]) 4673 else: d.append([]) 4674 # and replace 4675 self.data = d 4676 4677 # self.mdata 4678 lorig = self.mdata 4679 if lorig != None: 4680 # if length mis-match, skip (probably empty) 4681 if len(lorig) == nold: 4682 if self.verb > 3: print("++ padding mdata") 4683 d = [] 4684 # 1-based counting over runs 4685 for rind in rlist: 4686 # if > 0, get a full copy of an old run, else use an empty one 4687 if rind > 0: d.append(lorig[rind-1][:]) 4688 else: d.append([]) 4689 # and replace 4690 self.mdata = d 4691 4692 # self.alist 4693 lorig = self.alist 4694 if lorig != None: 4695 # if length mis-match, skip (probably empty) 4696 if len(lorig) == nold: 4697 if self.verb > 3: print("++ padding alist") 4698 d = [] 4699 # 1-based counting over runs 4700 for rind in rlist: 4701 # if > 0, get a full copy of an old run, else use an empty one 4702 if rind > 0: d.append(lorig[rind-1]) 4703 else: d.append(2) # 2 '*' for empty run 4704 # and replace 4705 self.alist = d 4706 4707 # self.run_lens 4708 lorig = self.run_lens 4709 if lorig != None: 4710 # if length mis-match, skip (probably empty) 4711 if len(lorig) == nold: 4712 if self.verb > 3: print("++ padding run_lens") 4713 d = [] 4714 # 1-based counting over runs 4715 for rind in rlist: 4716 # if > 0, get a full copy of an old run, else use an empty one 4717 if rind > 0: d.append(lorig[rind-1]) 4718 else: d.append(0) # empty run 4719 # and replace 4720 self.run_lens = d 4721 4722 if self.verb > 3: self.show("select_runs") 4723 4724 return 0 4725 4726 def show_duration_stats(self, per_run=0): 4727 """show min, mean, max, stdev for each column (unless col specified)""" 4728 4729 # if reasonable, show the file name 4730 fname = '' 4731 if self.verb: 4732 if self.fname: fname = "%s: " % self.fname 4733 form = "min = %g, mean = %g, max = %g, stdev = %g" 4734 else: 4735 form = "%g %g %g %g" 4736 4737 # apply newdata as end times, whine if negative 4738 dlist = [] 4739 for rind, run in enumerate(self.mdata): 4740 if per_run: 4741 mmms = form % UTIL.min_mean_max_stdev([event[2] for event in run]) 4742 print(('run %d : ' % rind) + mmms) 4743 else: 4744 dlist.extend([event[2] for event in run]) 4745 4746 # maybe finish with across-run display 4747 if not per_run: 4748 fstr = form % UTIL.min_mean_max_stdev(dlist) 4749 print('%s%s' % (fname, fstr)) 4750 4751 return 0 4752 4753 def show(self, mesg=''): 4754 print(self.make_show_str(mesg=mesg)) 4755 4756 def make_show_str(self, mesg=''): 4757 if mesg: mstr = '%s : ' % mesg 4758 else: mstr = '' 4759 4760 # report data elements via their top-level lengths 4761 ldata = -1 4762 if self.data != None: 4763 ldata = len(self.data) 4764 lmdata = -1 4765 if self.mdata != None: 4766 lmdata = len(self.mdata) 4767 lclines = -1 4768 if self.clines != None: 4769 lclines = len(self.clines) 4770 lalist = -1 4771 if self.alist != None: 4772 lalist = len(self.alist) 4773 4774 mstr = "--- %sAfniData element ---\n" \ 4775 " name (ready?) : %s (%d)\n" \ 4776 " fname : %s\n" \ 4777 " len(data, mdata, clines, alist) : %d, %d, %d, %d\n" \ 4778 " nrows, ncols, is_rect() : %d, %d, %d\n" \ 4779 " row_lens : %s\n" \ 4780 " binary, empty : %d, %d\n" \ 4781 " married, mtype : %d, %d\n" \ 4782 " minlen, maxlen, dur_len : %d, %d, %g\n" \ 4783 " tr, nruns : %g, %d\n" \ 4784 " run_lens : %s\n" \ 4785 " write_dm, cormat_ready : %d, %d\n" \ 4786 " verb : %d\n" \ 4787 % ( mstr, self.name, self.ready, self.fname, 4788 ldata, lmdata, lclines, lalist, 4789 self.nrows, self.ncols, self.is_rect(), self.row_lens, 4790 self.binary, self.empty, self.married, self.mtype, 4791 self.minlen, self.maxlen, self.dur_len, 4792 self.tr, self.nruns, self.run_lens, 4793 self.write_dm, self.cormat_ready, 4794 self.verb) 4795 4796 return mstr 4797 4798def show_multi_isi_stats(adlist, run_lens, tr, verb=0): 4799 from afnipy import lib_timing as LT 4800 4801 nad = len(adlist) 4802 if nad == 0: 4803 print('** show_multi_isi_stats: no elements to list') 4804 return 1 4805 4806 AD0 = adlist[0] 4807 MT = LT.AfniTiming(mdata=AD0.mdata) 4808 for ind in range(1, nad): 4809 tt = LT.AfniTiming(mdata=adlist[ind].mdata) 4810 MT.extend_rows(tt) 4811 4812 if verb: 4813 MT.show('multistim timing') 4814 4815 MT.show_isi_stats(mesg='%d elements'%nad, run_len=run_lens, tr=tr) 4816 4817if __name__ == '__main__': 4818 print('** this is not a main module') 4819 sys.exit(1) 4820 4821 4822