1#!/usr/local/bin/python3.8 2 3# python3 status: compatible 4 5# system libraries 6import sys, os 7 8# AFNI libraries 9from afnipy import option_list as OL 10from afnipy import afni_util as UTIL # not actually used, but probably will be 11from afnipy import lib_subjects as SUBJ 12 13# ---------------------------------------------------------------------- 14# globals 15 16g_help_string = """ 17============================================================================= 18gen_group_command.py - generate group analysis command scripts 19 20 purpose: ~1~ 21 22 Quickly generate group analysis command scripts by parsing wildcard-based 23 lists of input datasets. 24 25 1. generate group commands: 3dttest++, 3dMEMA, 3dANOVA2, 3dANOVA3 26 2. generate generic commands 27 3. todo (maybe): 3dttest, GroupAna (or maybe not) 28 29 This program is to assist in writing group commands. The hardest part (or 30 most tedious) is generally listing datasets and such, particularly including 31 sub-brick selection, and that is the main benefit of using this program. 32 33 If used without sufficient options (which might be typical), the generated 34 commands will not be complete (e.g. they might fail). So either provide 35 sufficient passed options via -options or plan to edit the resulting script. 36 37 If -write_script is not given, the command is written to stdout. 38 39 ** NOTE: this program expects one dataset per subject. Single condition 40 volumes are accessed using sub-brick selectors via -subs_betas 41 and possbily -subs_tstats. 42 43 This program can parse subject IDs from dataset names when the IDs are the 44 varying part of dataset names (e.g. stats_subj1234+tlrc.HEAD), as in: 45 46 gen_group_command.py -command 3dttest++ \\ 47 -dsets stats*+tlrc.HEAD 48 49 or when the subject IDs are the varying part of the directory names (while 50 the actual file names are identical), as in: 51 52 gen_group_command.py -command 3dttest++ \\ 53 -dsets subject_results/*/*.results/stats+tlrc.HEAD 54 55 56 Generic commands do not need to be part of AFNI. Perhaps one just wants 57 an orderly and indented list of file names to be part of a bigger script. 58 consider: 59 60 gen_group_command.py -command ls -dsets group_results/OL*D 61 62 or perhaps using 3dTcat to collect a sub-brick from each subject: 63 64 gen_group_command.py -command 3dTcat -subs_betas 'Arel#0_Coef' \\ 65 -dsets group_results/OL*D 66 67------------------------------------------ 68examples (by program) ~1~ 69 70 A. 3dttest++ (not 3dttest) ~2~ 71 72 Note: these commands apply to the sample group data under 73 AFNI_data6/group_results. 74 75 * Note: The 3dttest++ program defaults to setA minus setB, which is the 76 opposite of 3dttest and 3dMEMA (though it might be more natural). 77 The direction of the test can be further specified using either 78 -AminusB or -BminusA, which is always included in the resulting 79 command if there are 2 sets of data. 80 81 This program will always supply one of -AminusB or -BminusA, to be 82 clear. If the user does not provide one, -AminusB will be used. 83 84 Note also that 3dttest uses sub-brick labels which should make 85 this clear. 86 87 1. the most simple case, providing just the datasets ~3~ 88 89 The most simple case, providing just the datasets. The subject IDs 90 will be extracted from the dataset names. Since no sub-bricks are 91 provided, the betas will default to sub-brick 0 and the test will be 92 the mean compared with 0. 93 94 gen_group_command.py -command 3dttest++ \\ 95 -dsets REML*.HEAD 96 97 2. specifying set labels and beta weights for a 2-sample t-test ~3~ 98 99 Specify the sub-bricks and set labels to compare Vrel vs. Arel. 100 Write the command to the file cmd.tt++.2. 101 102 gen_group_command.py -command 3dttest++ \\ 103 -write_script cmd.tt++.2 \\ 104 -prefix tt++.2_V-A \\ 105 -dsets REML*.HEAD \\ 106 -set_labels Vrel Arel \\ 107 -subs_betas 'Vrel#0_Coef' 'Arel#0_Coef' 108 109 3. request a paired t-test and apply a mask ~3~ 110 111 gen_group_command.py -command 3dttest++ \\ 112 -write_script cmd.tt++.3 \\ 113 -prefix tt++.3_V-A_paired \\ 114 -dsets REML*.HEAD \\ 115 -set_labels Vrel Arel \\ 116 -subs_betas 'Vrel#0_Coef' 'Arel#0_Coef' \\ 117 -options \\ 118 -paired -mask mask+tlrc 119 120 4. include options specific to 3dttest++ (not gen_group_command.py) ~3~ 121 122 Exclude voxels that are identically zero across more than 20% of the 123 input datasets (presumably masked at the single subject level). 124 Convert output directly to z, since the DOF will vary across space. 125 126 gen_group_command.py -command 3dttest++ \\ 127 -write_script cmd.tt++.4 \\ 128 -prefix tt++.4_V-A_zskip \\ 129 -dsets REML*.HEAD \\ 130 -set_labels Vrel Arel \\ 131 -subs_betas 'Vrel#0_Coef' 'Arel#0_Coef' \\ 132 -options \\ 133 -zskip 0.8 -toz 134 135 5. including covariates and related options ~3~ 136 137 Use covariates to account for a sex difference. We might encode 138 females as 0 and males as 1 to get an intercept (main effect) that 139 applies to females (if we do not do any centering). However, we 140 want a main effect for the average between males and females, and 141 therefore have used -1 for males and +1 for females. Add NONE 142 for centering so that 3dttest++ does not do any. 143 144 Females have subject indices: 0, 1, 2, 3 and 5. 145 Males have subject indices: 4 and 6 through 9 (the last). 146 147 gen_group_command.py -command 3dttest++ \\ 148 -write_script cmd.tt++.5 \\ 149 -prefix tt++.5_covary \\ 150 -dsets data/OLSQ*.HEAD \\ 151 -subs_betas 'Vrel#0_Coef' \\ 152 -options \\ 153 -covariates sex_encode.txt \\ 154 -center NONE 155 156 157 6. specify index lists to restrict applied subject datasets ~3~ 158 159 Use -dset_index0_list to compare female subjects to males. 160 Both subject types are in the same directory (10 subjects total). 161 So the -dsets options will both specify the same list, which will 162 then be paired down via -dset_index0_list to indicate only females 163 and only males. 164 165 Females have subject indices: 0, 1, 2, 3 and 5. 166 Males have subject indices: 4 and 6 through 9 (the last). 167 168 gen_group_command.py -command 3dttest++ \\ 169 -write_script cmd.tt++.6 \\ 170 -prefix tt++.6_F-M \\ 171 -dsets data/OLSQ*.HEAD \\ 172 -dset_index0_list '0..3,5' \\ 173 -dsets data/OLSQ*.HEAD \\ 174 -dset_index0_list '4,6..$' \\ 175 -set_labels female male \\ 176 -subs_betas 'Vrel#0_Coef' 177 178 179 7. specify applied subjects via subject ID lists ~3~ 180 181 For BIDS, adjust subject IDs and get group lists from text files, 182 group1_subjects.txt and group2_subjects.txt. 183 184 gen_group_command.py \\ 185 -command 3dttest++ \\ 186 -write_script cmd.tt++.7 \\ 187 -prefix tt++.7_F-M \\ 188 -dsets sub-*/*.results/stats.sub*REML+tlrc.HEAD \\ 189 -dset_sid_list `cat group1_subjects.txt` \\ 190 -dsets sub-*/*.results/stats.sub*REML+tlrc.HEAD \\ 191 -dset_sid_list `cat group2_subjects.txt` \\ 192 -set_labels horses rabbits \\ 193 -subs_betas 'carrots#0_Coef' 194 195 196 See "3dttest++ -help" for details on its options. 197 198 -------------------- 199 200 B. 3dMEMA ~2~ 201 202 Note: these commands apply to the sample group data under 203 AFNI_data6/group_results. 204 205 Note: As with 3dttest, group comparisons are done as the second set minus 206 the first set. 207 208 209 1. most simple case, providing only datasets ~3~ 210 211 The most simple case, providing just the datasets. The subject IDs 212 will be extracted from the dataset names. Since no sub-bricks are 213 provided, the betas will be 0 and t-stats will be 1. 214 215 gen_group_command.py -command 3dMEMA \\ 216 -dsets REML*.HEAD 217 218 2. getting separate groups via directories ~3~ 219 220 This does not quite apply to AFNI_data6. Assuming there are 2 group 221 directories, write a 2-sample command. 222 223 gen_group_command.py -command 3dMEMA \\ 224 -write_script cmd.mema.2 \\ 225 -dsets groupA/REML*.HEAD \\ 226 -dsets groupB/REML*.HEAD 227 228 3. restrict subject datasets via an index list ~3~ 229 230 Run 3dMEMA, but restrict the subjects to partial lists from within 231 an entire list. This applies -dset_index0_list (or the sister 232 -dset_index1_list). 233 234 # assume these 9 subjects represent all under the 'data' dir 235 set subjects = ( AA BB CC DD EE FF GG HH II ) 236 237 a. Do a simple test on subjects AA, HH, II and FF. Indices are: 238 0-based: 0, 7, 8, 5 (AA=0, ..., II=8) 239 1-based: 1, 8, 9, 6 (AA=1, ..., II=9) 240 241 gen_group_command.py -command 3dMEMA \\ 242 -write_script cmd.mema.3a \\ 243 -dsets data/REML*.HEAD \\ 244 -dset_index0_list '0,7,8,5' 245 246 b. Do a test on sub-lists of subjects. 247 248 gen_group_command.py -command 3dMEMA \\ 249 -write_script cmd.mema.3b \\ 250 -dsets data/REML*.HEAD \\ 251 -dset_index0_list '0,7,8,5' \\ 252 -dsets data/REML*.HEAD \\ 253 -dset_index0_list '3,4,6,9' \\ 254 -subs_betas 'Arel#0_Coef' \\ 255 -subs_tstats 'Arel#0_Tstat' 256 257 See "3dMEMA -help" for details on the extra options. 258 259 -------------------- 260 261 C. 3dANOVA2 ~2~ 262 263 Note: these commands apply to the sample group data under 264 AFNI_data6/group_results. 265 266 Note: it seems better to create the script without any contrasts, and 267 add them afterwards (so the user can format well). However, if 268 no contrasts are given, the program will add 1 trivial one. 269 270 271 1. basic example, with datasets and volume indices ~3~ 272 273 The most simple case, providing just the datasets and a list of 274 sub-bricks. 275 276 gen_group_command.py -command 3dANOVA2 \\ 277 -dsets OLSQ*.HEAD \\ 278 -subs_betas 0 1 279 280 2. get more useful: ~3~ 281 282 - apply with a directory 283 - specify a script name 284 - specify a dataset prefix for the 3dANOVA2 command 285 - use labels for sub-brick indices 286 - specify a simple contrast 287 288 gen_group_command.py -command 3dANOVA2 \\ 289 -write_script cmd.A2.2 \\ 290 -prefix outset.A2.2 \\ 291 -dsets AFNI_data6/group_results/REML*.HEAD \\ 292 -subs_betas 'Vrel#0_Coef' 'Arel#0_Coef' \\ 293 -options \\ 294 -adiff 1 2 VvsA 295 296 -------------------- 297 298 D. 3dANOVA3 ~2~ 299 300 Note: these commands apply to the sample group data under 301 AFNI_data6/group_results. 302 303 Note: it seems better to create the script without any contrasts, and 304 add them afterwards (so the user can format well). However, if 305 no contrasts are given, the program will add 2 trivial ones, 306 just for a starting point. 307 308 Note: this applies either -type 4 or -type 5 from 3dANOVA3. 309 See "3dANOVA3 -help" for details on the types. 310 311 The user does not specify type 4 or 5. 312 313 type 4: there should be one -dsets option and a -factors option 314 type 5: there should be two -dsets options and no -factor 315 316 1. 3dANOVA3 -type 4 : simple ~3~ 317 318 This is a simple example of a 2-way factorial ANOVA (color by image 319 type), across many subjects. The colors are pink and blue, while the 320 images are of houses, faces and donuts. So there are 6 stimulus types 321 in this 2 x 3 design: 322 323 pink house pink face pink donut 324 blue house blue face blue donut 325 326 Since those were the labels given to 3dDeconvolve, the beta weights 327 will have #0_Coef appended, as in pink_house#0_Coef. Note that in a 328 script, the '#' character will need to be quoted. 329 330 There is only one set of -dsets given, as there are no groups. 331 332 gen_group_command.py -command 3dANOVA3 \\ 333 -dsets OLSQ*.HEAD \\ 334 -subs_betas \\ 335 "pink_house#0_Coef" "pink_face#0_Coef" "pink_donut#0_Coef" \\ 336 "blue_house#0_Coef" "blue_face#0_Coef" "blue_donut#0_Coef" \\ 337 -factors 2 3 338 339 2. 3dANOVA3 -type 4 : more useful ~3~ 340 341 Get more useful: 342 - apply with an input data directory 343 - specify a script name 344 - specify a dataset prefix for the 3dANOVA3 command 345 - specify simple contrasts 346 347 gen_group_command.py -command 3dANOVA3 \\ 348 -write_script cmd.A3.2 \\ 349 -prefix outset.A3.2 \\ 350 -dsets AFNI_data6/group_results/OLSQ*.HEAD \\ 351 -subs_betas \\ 352 "pink_house#0_Coef" "pink_face#0_Coef" "pink_donut#0_Coef" \\ 353 "blue_house#0_Coef" "blue_face#0_Coef" "blue_donut#0_Coef" \\ 354 -factors 2 3 \\ 355 -options \\ 356 -adiff 1 2 pink_vs_blue \\ 357 -bcontr -0.5 -0.5 1.0 donut_vs_house_face 358 359 3. 3dANOVA3 -type 5 : simple, with 2 groups ~3~ 360 361 Here is a simple case, providing just 2 groups of datasets and a list 362 of sub-bricks. 363 364 gen_group_command.py -command 3dANOVA3 \\ 365 -dsets OLSQ*.HEAD \\ 366 -dsets REML*.HEAD \\ 367 -subs_betas 0 1 368 369 4. 3dANOVA3 -type 5 : more detailed ~3~ 370 371 Get more useful: 372 - apply with an input data directory 373 - specify a script name 374 - specify a dataset prefix for the 3dANOVA3 command 375 - use labels for sub-brick indices 376 - specify simple contrasts 377 378 gen_group_command.py -command 3dANOVA3 \\ 379 -write_script cmd.A3.4 \\ 380 -prefix outset.A3.2 \\ 381 -dsets AFNI_data6/group_results/OLSQ*.HEAD \\ 382 -dsets AFNI_data6/group_results/REML*.HEAD \\ 383 -subs_betas 'Vrel#0_Coef' 'Arel#0_Coef' \\ 384 -options \\ 385 -adiff 1 2 OvsR \\ 386 -bdiff 1 2 VvsA 387 388 -------------------- 389 390 E. generic/other programs ~2~ 391 392 These commands apply to basically any program, as specified. Options 393 may be provided, along with 1 or 2 sets of data. If provided, the 394 -subs_betas selectors will be applied. 395 396 This might be useful for simply making part of a longer script, where 397 the dataset names are explicit. 398 399 400 1. very simple demonstration, for just an 'ls' command ~3~ 401 402 Perhaps a fairly useless example with 'ls', just for demonstration. 403 404 gen_group_command.py -command ls -dsets group_results/OL*D 405 406 2. using 3dTcat to collect a sub-brick from each subject ~3~ 407 408 gen_group_command.py -command 3dTcat -subs_betas 'Arel#0_Coef' \\ 409 -dsets group_results/OL*D 410 411 3. including 2 sets of subjects, with a different sub-brick per set ~3~ 412 413 gen_group_command.py -command 3dTcat -subs_betas 0 1 \\ 414 -dsets group_results/OLSQ*D \\ 415 -dsets group_results/REML*D 416 417 4. 2 sets of subjects ~3~ 418 419 Datasets in different directories, and with different sub-brick 420 selectors, along with: 421 422 - a script name (to write the script to a text file) 423 - a -prefix 424 - options for the command (just 1 in this case) 425 - common sub-brick selectors for dataset lists 426 427 gen_group_command.py -command 3dMean \\ 428 -write_script cmd.3dmean.txt \\ 429 -prefix aud_vid_stdev \\ 430 -options -stdev \\ 431 -subs_betas 'Arel#0_Coef' \\ 432 -dsets group_results/OLSQ*D \\ 433 -dsets group_results/REML*D 434 435------------------------------------------ 436command-line options: ~1~ 437------------------------------------------ 438terminal options: ~2~ 439 440 -help : show this help 441 -hist : show module history 442 -show_valid_opts : list valid options 443 -ver : show current version 444 445required parameters: ~2~ 446 447 -command COMMAND_NAME : resulting command, such as 3dttest++ 448 449 The current list of group commands is: 3dttest++, 3dMEMA, 3dANOVA2, 450 3dANOVA3. 451 452 3dANOVA2: applied as -type 3 only (factor x subjects) 453 3dANOVA3: -type 4: condition x condition x subject 454 (see -factors option) 455 -type 5: group x condition x subject 456 457 -dsets datasets ... : list of datasets 458 459 Each use of this option essentially describes one group of subjects. 460 All volumes for a given subject should be in a single dataset. 461 462 This option can be used multiple times, once per group. 463 464other options: ~2~ 465 466 -dset_sid_list SID SID ... : restrict -dsets datasets to this SID list 467 468 In some cases it is easy to use a wildcard to specify all datasets via 469 -dsets, but where subject groups would not be partitioned that way. 470 For example, you have a list of subjects to apply, per group, but no 471 way to separate them with a wildcard (e.g. in a BIDS tree, with no 472 group directories). 473 474 Consider this example: 475 476 -subj_prefix sub- \\ 477 -dsets sub-*/*.results/stats.sub*REML+tlrc.HEAD \\ 478 -dset_sid_list sub-0* 479 480 or make 2 subject lists, each starting with all subjects, but with 481 group lists contained in text files: 482 483 -subj_prefix sub- \\ 484 -dsets sub-*/*.results/stats.sub*REML+tlrc.HEAD \\ 485 -dset_sid_list `cat group1_subjects.txt` \\ 486 -dsets sub-*/*.results/stats.sub*REML+tlrc.HEAD \\ 487 -dset_sid_list `cat group2_subjects.txt` \\ 488 489 -dset_index0_list values... : restrict -dsets datasets to this 0-based list 490 -dset_index1_list values... : restrict -dsets datasets to this 1-based list 491 492 In some cases it is easy to use a wildcard to specify datasets via 493 -dsets, but there may be a grouping of subjects within that list. 494 For example, if both males and females are in the list of datasets 495 provided by -dsets, and if one wants a comparison between those 2 496 groups, then a pair of -dset_index0_list could be specified (1 for 497 each -dset) option to list which are the females and males. 498 499 Consider this example: 500 501 -dsets all/stats.*.HEAD \\ 502 -dset_index0_list '0..5,10..15' \\ 503 -dsets all/stats.*.HEAD \\ 504 -dset_index0_list '6..9,16..$' \\ 505 506 Note that -dsets is used twice, with IDENTICAL lists of datasets. 507 The respective -dset_index0_list options then restrict those lists to 508 0-based index lists, one for females, the other for males. 509 510 * One must be careful to get the indices correct, so check the output 511 command script to be sure the correct subjects are in each group. 512 513 The difference between -dset_index0_list and -dset_index1_list is just 514 that the former is a 0-based list (such as is used by AFNI programs), 515 while the latter is 1-based (such as is used by tcsh). A 0-based list 516 begins counting at 0 (as in offsets), while a list 1-based starts at 1. 517 Since use of either makes sense, both are provided. 518 519 For example, these options are equivalent: 520 521 -dset_index0_list 0,5..8 522 -dset_index1_list 1,6..9 523 524 The format for these index lists is the same as for AFNI sub-brick 525 selection. 526 527 -factors NF1 NF2 ... : list of factor levels, per condition 528 529 example: -factors 2 3 530 531 This option is currently only for '3dANOVA3 -type 4', which is a 532 condition x condition x subject test. It is meant to parse the 533 -subs_betas option, which lists all sub-bricks input to the ANOVA. 534 535 Assuming condition A has nA levels, and B has nB (2 and 3 in the 536 above example), then this option (applied '-factors nA nB', and 537 -subs_betas) would take nA * nB parameters (for the cross product of 538 factor A and factor B levels). 539 The betas should be specified in A major order, as in: 540 541 -subs_betas A1B1_name A1B2_name ... A1BnB A2B1 A2B2 ... AnABnB_name 542 543 or as in the 2 x 3 case: 544 545 -subs_betas A1B1 A1B2 A1B3 A2B1 A2B2 A2B3 -factors 2 3 546 547 e.g. for pink/blue x house/face/donut, output be 3dDeconvolve 548 (i.e. each betas probably has #0_Coef attached) 549 550 -subs_betas \\ 551 "pink_house#0_Coef" "pink_face#0_Coef" "pink_donut#0_Coef" \\ 552 "blue_house#0_Coef" "blue_face#0_Coef" "blue_donut#0_Coef" \\ 553 -factors 2 3 \\ 554 555 Again, these factor combination names should be either sub-brick labels 556 or indices (labels are suggested, to avoid confusion). 557 558 See the example with '3dANOVA3 -type 4' as part of example D, above. 559 See also -subs_betas. 560 561 -keep_dirent_pre : keep directory entry prefix 562 563 Akin to -subj_prefix, this flag expands the subject prefix list to 564 include everything up to the beginning of the directory names (at 565 the level that varies across input datasets). 566 567 Example 1: 568 datasets: 569 subj.FP/betas+tlrc subj.FR/betas+tlrc subj.FT/betas+tlrc 570 subj.FV/betas+tlrc subj.FW/betas+tlrc subj.FX/betas+tlrc 571 subj.FY/betas+tlrc subj.FZ/betas+tlrc 572 573 The default subject IDs would be: 574 P R T V W X Y Z 575 576 When using -keep_dirent_pre, subject IDs would be: 577 subj.FP subj.FR subj.FT subj.FV subj.FW subj.FX subj.FY subj.FZ 578 579 Note that these IDs come at the directory level, since the dataset 580 names do not vary. 581 582 Example 2: 583 datasets: 584 subj.FP/OLSQ.FP.betas+tlrc subj.FR/OLSQ.FR.betas+tlrc 585 subj.FT/OLSQ.FT.betas+tlrc subj.FV/OLSQ.FV.betas+tlrc 586 subj.FW/OLSQ.FW.betas+tlrc subj.FX/OLSQ.FX.betas+tlrc 587 subj.FY/OLSQ.FY.betas+tlrc subj.FZ/OLSQ.FZ.betas+tlrc 588 589 The default subject IDs would be: 590 P R T V W X Y Z 591 592 When using -keep_dirent_pre, subject IDs would be: 593 OLSQ.FP OLSQ.FR OLSQ.FT OLSQ.FV OLSQ.FW OLSQ.FX OLSQ.FY OLSQ.FZ 594 595 Note that these IDs come at the dataset level, since the dataset 596 names vary. 597 598 -hpad PAD : pad subject prefix by PAD chars toward header 599 600 Akin to -subj_prefix and -tpad, this flag expands the subject prefix 601 list to include PAD extra characters toward the beginning. 602 603 See also -tpad. 604 605 -tpad PAD : pad subject prefix by PAD chars toward tail 606 607 Akin to -subj_prefix and -hpad, this flag expands the subject prefix 608 list to include PAD extra characters toward the beginning. 609 610 See also -hpad. 611 612 -options OPT1 OPT2 ... : list of options to pass along to result 613 614 The given options will be passed directly to the resulting command. If 615 the -command is 3dMEMA, say, these should be 3dMEMA options. This 616 program will not evaluate or inspect the options, but will put them at 617 the end of the command. 618 619 -prefix PREFIX : apply as COMMAND -prefix 620 -set_labels LAB1 LAB2 ... : labels corresponding to -dsets entries 621 -subj_prefix PREFIX : prefix for subject names (3dMEMA) 622 -subj_suffix SUFFIX : suffix for subject names (3dMEMA) 623 -subs_betas B0 B1 : sub-bricks for beta weights (or similar) 624 625 If this option is not given, sub-brick 0 will be used. The entries 626 can be either numbers or labels (which should match what is seen in 627 the afni GUI, for example). 628 629 If there are 2 -set_labels, there should be 2 betas (or no option). 630 631 -subs_tstats T0 T1 : sub-bricks for t-stats (3dMEMA) 632 633 If this option is not given, sub-brick 1 will be used. The entries can 634 be either numbers or labels (which should match what is seen in the 635 afni GUI, for example). 636 637 This option applies only to 3dMEMA currently, and in that case, its use 638 should match that of -subs_betas. 639 640 See also -subs_betas. 641 642 -type TEST_TYPE : specify the type of test to perform 643 644 The test type may depend on the given command, but generally implies 645 there are multiple sets of values to compare. Currently valid tests 646 are (for the given program): 647 648 3dMEMA: paired, unpaired 649 650 If this option is not applied, a useful default will be chosen. 651 652 -verb LEVEL : set the verbosity level 653 654 -write_script FILE_NAME : write command script to FILE_NAME 655 656 If this option is given, the command will be written to the specified 657 file name. Otherwise, it will be written to the terminal window. 658 659----------------------------------------------------------------------------- 660R Reynolds October 2010 661============================================================================= 662""" 663 664g_history = """ 665 gen_group_command.py history: 666 667 0.0 Sep 09, 2010 - initial version 668 0.1 Oct 25, 2010 - handle some 3dMEMA cases 669 0.2 Oct 26, 2010 - MEMA updates 670 0.3 Nov 08, 2010 - can generate 3dttest++ commands 671 0.4 Jun 15, 2011 - if constant dset names, extract SIDs from dir names 672 (done for R Momenan) 673 0.5 Jun 27, 2011 674 - added -dset_index0_list/-dset_index1_list options (for R Momenan) 675 - ttest++ and MEMA commands now apply directories to datasets 676 - changed Subject.atrs to be VarsObject instance, not dictionary 677 0.6 Jun 22, 2012 678 - added commands 3dANOVA2 and 3dANOVA3 679 - added -factors for 3dANOVA3 -type 4 680 0.7 Jun 25, 2012 - added help for -factors and 3dANOVA3 -type 4 examples 681 0.8 Sep 04, 2012 - fixed error message 682 0.9 Oct 03, 2012 - some options do not allow dashed parameters 683 0.10 Oct 30, 2013 - added -keep_dirent_pre 684 0.11 Apr 24, 2015 - tiny help update (example) 685 0.12 Aug 12, 2015 - allow generic (unknown) commands (via -command) 686 0.13 Mar 29, 2016 - 3dMEMA now requires paired test to be via input contrast 687 1.0 Dec 27, 2017 - python3 compatible 688 1.1 Feb 26, 2019 689 - added -dset_sid_list to specify subject list, like -dset_index0_list 690 - added -hpad and -tpad opts 691 - use less indentation to tighten 3dttest++ command (do others?) 692 1.2 Mar 5, 2019 693 - show subject counts 694 - change max line len and whether data dir vars are used 695 - no require on restricted subjects 696 1.3 Jul 30, 2019 - sphinx help update 697""" 698 699g_version = "gen_group_command.py version 1.3 July 30, 2019" 700 701 702class CmdInterface: 703 """interface class for getting commands from SubjectList class 704 """ 705 706 def __init__(self, verb=1): 707 # main variables 708 self.status = 0 # exit value 709 self.valid_opts = None 710 self.user_opts = None 711 712 # general variables 713 self.command = '' # program name to make command for 714 self.ttype = None # test type (e.g. paired) 715 self.comp_dir = '-AminusB' # contrast direction (or -BminusA) 716 self.prefix = None # prefix for command result 717 self.write_script = None # file to write output to (else stdout) 718 self.betasubs = None # list of beta weight sub-brick indices 719 self.tstatsubs = None # list of t-stat sub-brick indices 720 self.lablist = None # list of set labels 721 self.factors = [] # list of factors of each type 722 self.hpad = 0 # hpad for list_minus_glob_form 723 self.tpad = 0 # tpad for list_minus_glob_form 724 725 self.subj_prefix = '' # prefix for each subject ID 726 self.subj_suffix = '' # suffix for each subject ID 727 self.dent_pre = 0 # flag: keep dir entry prefix 728 self.verb = verb 729 730 # lists 731 self.options = [] # other command options 732 self.slist = [] # list of SubjectList elements 733 self.dsets = [] # list of lists of filenames 734 self.index0_list = [] # 0-based sub-list of 'dsets' 735 self.index1_list = [] # 1-based sub-list of 'dsets' 736 self.sid_apply = [] # lists of subject IDs to apply 737 self.sid_omit = [] # lists of subject IDs to omit 738 739 # initialize valid_opts 740 self.init_options() 741 742 def show(self): 743 print("---------------------------- setup -----------------------------") 744 print("command : %s" % self.command) 745 print("test type : %s" % self.ttype) 746 print("prefix : %s" % self.prefix) 747 print("write_script : %s" % self.write_script) 748 print("beta sub-bricks : %s" % self.betasubs) 749 print("tstat sub-bricks : %s" % self.tstatsubs) 750 print("label list : %s" % self.lablist) 751 print("subject prefix : %s" % self.subj_prefix) 752 print("subject suffix : %s" % self.subj_suffix) 753 print("dirent prefix : %s" % self.dent_pre) 754 print("verb : %s" % self.verb) 755 756 print("options : %s" % self.options) 757 print("subject list(s) : %s" % self.slist) 758 print("index0_list : %s" % self.index0_list) 759 print("index1_list : %s" % self.index1_list) 760 print("datasets : %s" % self.dsets) # last 761 762 if self.verb > 3: 763 print("status : %s" % self.status) 764 self.user_opts.show(mesg="user options : ") 765 print("----------------------------------------------------------------") 766 767 def init_options(self): 768 self.valid_opts = OL.OptionList('valid opts') 769 770 # terminal options 771 self.valid_opts.add_opt('-help', 0, [], \ 772 helpstr='display program help') 773 self.valid_opts.add_opt('-hist', 0, [], \ 774 helpstr='display the modification history') 775 self.valid_opts.add_opt('-show_valid_opts', 0, [],\ 776 helpstr='display all valid options') 777 self.valid_opts.add_opt('-ver', 0, [], \ 778 helpstr='display the current version number') 779 780 # required parameters 781 self.valid_opts.add_opt('-command', 1, [], 782 helpstr='specify the program used in the output command') 783 self.valid_opts.add_opt('-dsets', -1, [], okdash=0, 784 helpstr='specify a list of input datasets') 785 786 # other options 787 self.valid_opts.add_opt('-AminusB', 0, [], 788 helpstr='apply 3dttest++ test as set A minus set B') 789 self.valid_opts.add_opt('-BminusA', 0, [], 790 helpstr='apply 3dttest++ test as set B minus set A') 791 self.valid_opts.add_opt('-dset_index0_list', -1, [], okdash=0, 792 helpstr='restrict dsets to 0-based index list') 793 self.valid_opts.add_opt('-dset_index1_list', -1, [], okdash=0, 794 helpstr='restrict dsets to 1-based index list') 795 self.valid_opts.add_opt('-dset_sid_list', -1, [], okdash=0, 796 helpstr='restrict dsets to these subject IDs') 797 self.valid_opts.add_opt('-dset_sid_omit_list', -1, [], okdash=0, 798 helpstr='remove these subject IDs from dsets') 799 self.valid_opts.add_opt('-factors', -1, [], okdash=0, 800 helpstr='num factors, per condition (probably 2 ints)') 801 self.valid_opts.add_opt('-hpad', 1, [], okdash=0, 802 helpstr='pad header by this length in subj IDs') 803 self.valid_opts.add_opt('-tpad', 1, [], okdash=0, 804 helpstr='pad tail by this length in subj IDs') 805 self.valid_opts.add_opt('-keep_dirent_pre', 0, [], 806 helpstr='keep directory entry prefix') 807 self.valid_opts.add_opt('-options', -1, [], 808 helpstr='specify options to pass to the command') 809 self.valid_opts.add_opt('-prefix', 1, [], 810 helpstr='specify output prefix for the command') 811 self.valid_opts.add_opt('-set_labels', -1, [], okdash=0, 812 helpstr='list of labels for each set of subjects') 813 self.valid_opts.add_opt('-subj_prefix', 1, [], 814 helpstr='specify prefix for each subject ID') 815 self.valid_opts.add_opt('-subj_suffix', 1, [], 816 helpstr='specify suffix for each subject ID') 817 self.valid_opts.add_opt('-subs_betas', -1, [], okdash=0, 818 helpstr='beta weight sub-bricks, one per subject list') 819 self.valid_opts.add_opt('-subs_tstats', -1, [], okdash=0, 820 helpstr='t-stat sub-bricks, one per subject list') 821 self.valid_opts.add_opt('-type', 1, [], 822 helpstr='specify the test type (e.g. paired)') 823 self.valid_opts.add_opt('-verb', 1, [], 824 helpstr='set the verbose level (default is 1)') 825 self.valid_opts.add_opt('-write_script', 1, [], 826 helpstr='specify file to write command into') 827 828 return 0 829 830 def process_options(self, argv=sys.argv): 831 832 # process any optlist_ options 833 self.valid_opts.check_special_opts(argv) 834 835 # process terminal options without the option_list interface 836 # (so that errors are not reported) 837 838 # if no arguments are given, apply -help 839 if len(argv) <= 1 or '-help' in argv: 840 print(g_help_string) 841 return 0 842 843 if '-hist' in argv: 844 print(g_history) 845 return 0 846 847 if '-show_valid_opts' in argv: 848 self.valid_opts.show('', 1) 849 return 0 850 851 if '-ver' in argv: 852 print(g_version) 853 return 0 854 855 # ============================================================ 856 # read options specified by the user 857 self.user_opts = OL.read_options(argv, self.valid_opts) 858 uopts = self.user_opts # convenience variable 859 if not uopts: return 1 # error condition 860 861 # ------------------------------------------------------------ 862 # require a list of files, at least 863 864 # ------------------------------------------------------------ 865 # process options, go after -verb first 866 867 val, err = uopts.get_type_opt(int, '-verb') 868 if val != None and not err: self.verb = val 869 870 for opt in uopts.olist: 871 872 # main options 873 if opt.name == '-AminusB': 874 self.comp_dir = opt.name 875 continue 876 877 if opt.name == '-BminusA': 878 self.comp_dir = opt.name 879 continue 880 881 if opt.name == '-command': 882 val, err = uopts.get_string_opt('', opt=opt) 883 if val == None or err: return 1 884 self.command = val 885 continue 886 887 if opt.name == '-dsets': 888 val, err = uopts.get_string_list('', opt=opt) 889 if val == None or err: return 1 890 self.dsets.append(val) # allow multiple such options 891 continue 892 893 if opt.name == '-dset_index0_list': 894 val, err = uopts.get_string_list('', opt=opt) 895 if val == None or err: return 1 896 self.index0_list.append(val) # allow multiple such options 897 continue 898 899 if opt.name == '-dset_index1_list': 900 val, err = uopts.get_string_list('', opt=opt) 901 if val == None or err: return 1 902 self.index1_list.append(val) # allow multiple such options 903 continue 904 905 if opt.name == '-dset_sid_list': 906 val, err = uopts.get_string_list('', opt=opt) 907 if val == None or err: return 1 908 self.sid_apply.append(val) # allow multiple such options 909 continue 910 911 if opt.name == '-dset_sid_omit_list': 912 val, err = uopts.get_string_list('', opt=opt) 913 if val == None or err: return 1 914 self.sid_omit.append(val) # allow multiple such options 915 continue 916 917 if opt.name == '-factors': 918 val, err = uopts.get_type_list(int, '', opt=opt) 919 if val == None or err: return 1 920 self.factors = val 921 continue 922 923 if opt.name == '-hpad': 924 val, err = uopts.get_type_opt(int, opt=opt) 925 if val == None or err: return 1 926 self.hpad = val 927 continue 928 929 if opt.name == '-tpad': 930 val, err = uopts.get_type_opt(int, opt=opt) 931 if val == None or err: return 1 932 self.tpad = val 933 continue 934 935 if opt.name == '-keep_dirent_pre': 936 self.dent_pre = 1 937 continue 938 939 if opt.name == '-options': 940 val, err = uopts.get_string_list('', opt=opt) 941 if val == None or err: return 1 942 self.options = val 943 continue 944 945 if opt.name == '-prefix': 946 val, err = uopts.get_string_opt('', opt=opt) 947 if val == None or err: return 1 948 self.prefix = val 949 continue 950 951 if opt.name == '-set_labels': 952 val, err = uopts.get_string_list('', opt=opt) 953 if val == None or err: return 1 954 self.lablist = val 955 continue 956 957 if opt.name == '-subj_prefix': 958 val, err = uopts.get_string_opt('', opt=opt) 959 if val == None or err: return 1 960 self.subj_prefix = val 961 continue 962 963 if opt.name == '-subj_suffix': 964 val, err = uopts.get_string_opt('', opt=opt) 965 if val == None or err: return 1 966 self.subj_suffix = val 967 continue 968 969 if opt.name == '-subs_betas': 970 val, err = uopts.get_string_list('', opt=opt) 971 if val == None or err: return 1 972 self.betasubs = val 973 continue 974 975 if opt.name == '-subs_tstats': 976 val, err = uopts.get_string_list('', opt=opt) 977 if val == None or err: return 1 978 self.tstatsubs = val 979 continue 980 981 if opt.name == '-type': 982 val, err = uopts.get_string_opt('', opt=opt) 983 if val == None or err: return 1 984 self.ttype = val 985 continue 986 987 if opt.name == '-verb': continue # already handled 988 989 if opt.name == '-write_script': 990 val, err = uopts.get_string_opt('', opt=opt) 991 if val == None or err: return 1 992 self.write_script = val 993 continue 994 995 # general options 996 997 # an unhandled option 998 print('** option %s not yet supported' % opt.name) 999 return 1 1000 1001 if self.verb > 2: self.show() 1002 1003 # process -dset_index_list_0 and _1 1004 if self.update_dset_lists(): return 1 1005 1006 return None 1007 1008 def update_dset_lists(self): 1009 """process and -dset_index0_list or -dset_index1_list options 1010 (do not allow both) 1011 """ 1012 1013 oname0 = '-dset_index0_list' 1014 oname1 = '-dset_index1_list' 1015 1016 uopts = self.user_opts 1017 1018 # if no list selectors, there is nothing to do 1019 if len(self.index0_list) == 0 and len(self.index1_list) == 0: return 0 1020 1021 # both option types is an error 1022 if len(self.index0_list) > 0 and len(self.index1_list) > 0: 1023 print('** cannot use both %s and %s' % (oname0, oname1)) 1024 return 1 1025 1026 otype = 0 1027 oname = oname0 1028 olist = self.index0_list 1029 nopt = len(olist) 1030 1031 if nopt == 0: 1032 otype = 1 1033 oname = oname1 1034 olist = self.index1_list 1035 nopt = len(olist) 1036 1037 # require one -dset_index option per -dsets option 1038 if nopt != len(self.dsets): 1039 print('** num -dset_indexX_list opts must match num -dsets opts' \ 1040 ' (%d != %d)' % (nopt, len(self.dsets))) 1041 return 1 1042 1043 new_dsets = [] 1044 for dind, dlist in enumerate(self.dsets): 1045 status, newlist = UTIL.restrict_by_index_lists(dlist, olist[dind], 1046 otype, nonempty=1, verb=self.verb) 1047 if status: 1048 print('** bad use of %s' % oname) 1049 return 1 1050 new_dsets.append(newlist) 1051 1052 self.dsets = new_dsets 1053 1054 def execute(self): 1055 1056 if not self.ready_for_action(): return 1 1057 1058 if self.verb > 1: 1059 print('-- make %s command with %d set(s) of dsets of length(s): %s' \ 1060 % (self.command, len(self.dsets), 1061 ', '.join([str(len(dlist)) for dlist in self.dsets]) )) 1062 1063 n_sid_apply = len(self.sid_apply) 1064 if n_sid_apply > 0 and not (n_sid_apply == len(self.dsets)): 1065 print("** num -dset_sid_list opts should match -dsets") 1066 return 1 1067 1068 n_sid_omit = len(self.sid_omit) 1069 if n_sid_omit > 0 and not (n_sid_omit == len(self.dsets)): 1070 print("** num -dset_sid_omit_list opts should match -dsets") 1071 return 1 1072 1073 # array of [total, post-restricted, post-removed], per list 1074 subj_count = [] 1075 # might deal with subject IDs and attributes later 1076 for ind, dlist in enumerate(self.dsets): 1077 slist = SUBJ.SubjectList(dset_l=dlist, verb=self.verb) 1078 if slist.status: return 1 1079 if slist.set_ids_from_dsets(prefix=self.subj_prefix, 1080 suffix=self.subj_suffix, 1081 hpad=self.hpad, 1082 tpad=self.tpad, 1083 dpre=self.dent_pre): 1084 print('** cannot set subject IDs from datasets') 1085 return 1 1086 scount = [len(slist.subjects)] # total 1087 1088 # possibly restrict subject lists to those chosen 1089 if n_sid_apply > 0: 1090 if slist.restrict_ids_to_dsets(self.sid_apply[ind], require=0): 1091 return 1 1092 scount.append(len(slist.subjects)) # restricted 1093 1094 # and possibly remove the undesirables 1095 if n_sid_omit > 0: 1096 if slist.remove_ids_from_dsets(self.sid_omit[ind], require=0): 1097 return 1 1098 scount.append(len(slist.subjects)) # restricted 1099 subj_count.append(scount) # and append current counts 1100 1101 # and store the list 1102 self.slist.append(slist) 1103 if self.verb > 2: slist.show("slist %d" % ind) 1104 1105 if self.verb > 1: 1106 print("subject counts:") 1107 print(" %-16s %-16s %-16s %-16s" \ 1108 % ('label', 'init nsubj', 'after restrict', 'after omit')) 1109 for scind, sc in enumerate(subj_count): 1110 if self.lablist and (len(self.lablist) == len(subj_count)): 1111 slab = self.lablist[scind] 1112 else: 1113 slab = 'slist_%d' % scind 1114 print(" %-16s %-16s %-16s %-16s" % (slab, sc[0], sc[1], sc[2])) 1115 print("") 1116 1117 cmd = None 1118 if self.command == '3dMEMA': 1119 cmd = self.get_mema_command() 1120 elif self.command == '3dttest++': 1121 cmd = self.get_ttpp_command() 1122 elif self.command == '3dANOVA2': 1123 cmd = self.get_anova2_command() 1124 elif self.command == '3dANOVA3': 1125 cmd = self.get_anova3_command() 1126 elif self.command: 1127 cmd = self.get_generic_command() 1128 else: 1129 print('** command not implemented: %s' % self.command) 1130 1131 # bail on failure, else wrap command 1132 if cmd == None: 1133 print('** failed making %s command' % self.command) 1134 return 1 1135 cmd = UTIL.add_line_wrappers(cmd, maxlen=100) 1136 1137 # either write to file or print 1138 if self.write_script: 1139 ofile = self.write_script 1140 if UTIL.write_text_to_file(ofile, cmd): 1141 print("** failed to write command to file '%s'" % ofile) 1142 return 1 1143 if self.verb > 0 and ofile != '-' and ofile != 'stdout': 1144 print('++ command written to file %s' % ofile) 1145 else: print(cmd) 1146 1147 def get_mema_command(self): 1148 if len(self.slist) > 1: s2 = self.slist[1] 1149 else: s2 = None 1150 if (self.betasubs != None and self.tstatsubs == None) or \ 1151 (self.betasubs == None and self.tstatsubs != None): 1152 print('** MEMA: -subs_betas and -subs_tstats must be used together') 1153 return None 1154 return self.slist[0].make_mema_command(set_labs=self.lablist, 1155 bsubs=self.betasubs, tsubs=self.tstatsubs, subjlist2=s2, 1156 prefix=self.prefix, ttype=self.ttype, options=self.options) 1157 1158 def get_generic_command(self): 1159 if len(self.slist) > 1: s2 = self.slist[1] 1160 else: s2 = None 1161 return self.slist[0].make_generic_command(self.command, 1162 bsubs=self.betasubs, subjlist2=s2, prefix=self.prefix, 1163 options=self.options) 1164 1165 def get_ttpp_command(self): 1166 if len(self.slist) > 1: s2 = self.slist[1] 1167 else: s2 = None 1168 return self.slist[0].make_ttestpp_command(set_labs=self.lablist, 1169 bsubs=self.betasubs, subjlist2=s2, prefix=self.prefix, 1170 comp_dir=self.comp_dir, options=self.options) 1171 1172 def get_anova2_command(self): 1173 """generate 3dANOVA2 command 1174 type 2: requires one group and one list of betas 1175 """ 1176 if self.betasubs == None: 1177 print('** missing required -subs_betas option for sub-brick list') 1178 return None 1179 1180 return self.slist[0].make_anova2_command( bsubs=self.betasubs, 1181 prefix=self.prefix, options=self.options, verb=self.verb) 1182 1183 def get_anova3_command(self): 1184 """generate 3dANOVA3 command 1185 type 5: requires 2 groups and one list of betas 1186 """ 1187 1188 return self.slist[0].make_anova3_command( bsubs=self.betasubs, 1189 prefix=self.prefix, subjlists=self.slist, options=self.options, 1190 factors=self.factors, verb=self.verb) 1191 1192 def help_mema_command(self): 1193 helpstr = """ 1194 3dMEMA command help: 1195 1196 This is for help in deciding which MEMA command format to use, which 1197 command parameters are needed, and how a command would be organized. 1198 1199 As with 3dttest, there are 3 basic ways to run 3dMEMA. 1200 1201 1. as a one-sample test (see example 1 from '3dMEMA -help') 1202 2. as a two-sample test (see example 3 from '3dMEMA -help') 1203 3. as a paired test (see example 2 from '3dMEMA -help') 1204 1205 1. For the one-sample test, the required inputs are the datasets. 1206 It is best to also supply corresponding beta and t-stat sub-brick 1207 indexes or labels as well. 1208 1209 minimum: 1210 1211 -dsets stats.*+tlrc.HEAD 1212 1213 suggested: 1214 1215 -subs_betas 0 -subs_tstats 1 1216 OR 1217 -subs_betas 'Vrel#0_Coef' -subs_tstats 'Vrel#0_Tstat' 1218 1219 2. 1220 """ 1221 1222 def help_datasets(self): 1223 helpstr = """ 1224 Dataset configuration and naming: 1225 When using this command generator, each subject should have all 1226 data (betas and t-stats) in a single dataset. Dataset names 1227 should preferably be consistent, varying over only subject ID 1228 codes (this is not a requirement, but makes life easier). 1229 1230 For example, this allows one to specify datasets at once (via a 1231 wildcard) and let the program sort it out. Consider using: 1232 1233 -dsets stats.*+tlrc.HEAD 1234 1235 This will expand alphabetically, with subject IDs coming from 1236 the part of the dataset names that varies. 1237 1238 """ 1239 1240 def ready_for_action(self): 1241 1242 ready = 1 1243 1244 if self.command == '': 1245 if self.verb > 0: print('** missing execution command') 1246 ready = 0 1247 1248 if len(self.dsets) < 1: 1249 if self.verb > 0: print('** missing datasets for command') 1250 ready = 0 1251 1252 return ready 1253 1254 def init_from_file(self, fname): 1255 """load a 1D file, and init the main class elements""" 1256 1257 self.status = 1 # init to failure 1258 adata = L1D.Afni1D(fname, name='1d_tool data', verb=self.verb) 1259 if not adata.ready: 1260 print("** failed to read 1D data from '%s'" % fname) 1261 return 1 1262 1263 if self.verb > 1: print("++ read 1D data from file '%s'" % fname) 1264 1265 self.ad = adata 1266 self.status = 0 1267 1268 return 0 1269 1270 def test(self, verb=3): 1271 print('------------------------ initial tests -----------------------') 1272 self.verb = verb 1273 # first try AFNI_data4, then regression data 1274 1275 print('------------------------ reset files -----------------------') 1276 1277 print('------------------------ should fail -----------------------') 1278 1279 print('------------------------ more tests ------------------------') 1280 1281 return None 1282 1283def main(): 1284 me = CmdInterface() 1285 if not me: return 1 1286 1287 rv = me.process_options() 1288 if rv != None: return rv 1289 1290 rv = me.execute() 1291 if rv != None: return rv 1292 1293 return me.status 1294 1295if __name__ == '__main__': 1296 sys.exit(main()) 1297 1298 1299