-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathoptimizeXGUI.m
executable file
·3189 lines (2751 loc) · 110 KB
/
optimizeXGUI.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
function optimizeXGUI(params)
% OPTIMIZEXGUI
%
% USAGE: optimizeXGUI([params])
%
% Run this function at the command line to open the GUI. Once the design
% search has ended, a time stamped folder will be created in the software
% directory. In the folder, details regarding the most efficient designs
% will be saved. This includes a .mat file with all of the design
% information, as well as separate txt/csv files for the most efficient
% designs. The text files will have 5 columns:
%
% column 1 -- trial #
% column 2 -- condition
% column 3 -- onsets (in secs)
% column 4 -- trial duration (secs)
% column 5 -- post-trial duration (secs)
%
%
% The GUI was created using the MATLAB File Exchange contribution
% "settingsdlg.m" (http://goo.gl/DFvcQ5), which is included as a
% subfunction. Some of the code is lifted or adapted from software
% developed by Russ Poldrack or Tor Wager.
% Copyright 2014 Bob Spunt, Pasadena, California
% California Institute of Technology, USA
% Contact: bobspunt@gmail.com
% Current Version: 20160104
%
if nargin==0
params = getParams;
elseif ~isstruct(params)
if iscell(params), params = char(params); end
try
load(params)
catch err
error_handler(err);
end
end
% =========================================================================
% Run defineX.m
% =========================================================================
params = defineX(params);
% =========================================================================
% Run optimizeX.m
% =========================================================================
try
optimizeX(params);
delete(findall(0, 'Name','optimizeXGUI Progress'));
catch err
error_handler(err);
delete(findall(0, 'Name','optimizeXGUI Progress'));
end
end
% =========================================================================
% *
% * SUBFUNCTIONS (PRIMARY)
% *
% =========================================================================
function params = getParams
% =========================================================================
% Request Major Settings from User
% =========================================================================
try
% cbintervalstr = {'None', '1st/2nd Half (2)', 'Tertiles (3)', 'Quartiles (4)', 'Quintiles (5)' , 'Sextile (6)', 'Septile (7)', 'Octile (8)', 'Decile (10)', 'Hexadecile (16)'};
% cbintervalnum = [1 2 3 4 5 6 7 8 10 16];
% cbintervalstr = {'None', '1st/2nd Half (2)', 'Tertiles (3) [possibly slow]', 'Quartiles (4) [possibly very slow]'};
% cbintervalnum = [1 2 3 4];
% {'Interval Counterbalancing'; 'counterBalanceInterval'}, cbintervalstr, ...
distOpt = {
'Rayleigh'
'Chisquare'
'Exponential'
};
[params, button] = settingsdlg(...
'title' , 'OptimizeX Settings', ...
'separator' , 'General Settings', ...
{'TR (s)'; 'TR'}, 2, ...
{'High-Pass Filter Cutoff (s)'; 'hpf'}, 100, ...'
'separator' , 'Task Settings', ...
{'N Conditions';'nconds'} , 4, ...
{'N Trials Per Condition';'trialsPerCond'} , '25 25 25 25', ...
{'Maximum Block Size'; 'maxRep'} , 3, ...
'separator' , 'Timing (s)', ...
{'Trial Duration'; 'trialDur'}, 2, ...
{'Mean ISI';'meanISI'} , 3, ...
{'Min ISI';'minISI'} , 2, ...
{'Max ISI';'maxISI'} , 6, ...
{'ISI Distribution'; 'distISI'} , distOpt, ...
{'Time before first trial'; 'restBegin'}, 10, ...
{'Time after last trial'; 'restEnd'}, 10, ...
'separator' , 'Optimization Settings', ...
{'N Designs to Save'; 'keep'}, 5, ...
{'N Generations to Run';'ngen'} , 100, ...
{'N Designs Per Generation';'gensize'} , 500, ...
{'Max Time to Run (minutes)';'maxtime'} , .5);
catch err
error_handler(err);
end
if strcmpi(button, 'cancel') || isempty(button), return; end % canceled
% params.counterBalanceInterval = cbintervalnum(strcmpi(cbintervalstr, params.counterBalanceInterval));
params.trialsPerCond = str2num(params.trialsPerCond);
params.counterBalanceInterval = 0;
if length(params.trialsPerCond)==1, params.trialsPerCond = repmat(params.trialsPerCond, 1, params.nconds);
elseif length(params.trialsPerCond)~=params.nconds
msg = sprintf('The number of entries in "N Trials Per Condition" does not match the number of conditions');
headsup(msg, 'Invalid Input', 1);
optimizeXGUI
end
if params.minISI > params.meanISI | params.maxISI < params.meanISI
msg = sprintf('The ISI values you''ve specified look odd: Min ISI cannot be greater than the Mean ISI, and the Max ISI cannot be less than the Mean ISI');
headsup(msg, 'Invalid Input', 1);
optimizeXGUI
end
% =========================================================================
% Contrast Specification and Weighting
% =========================================================================
[condata, button] = settingsdlg(...
'title' , 'Settings', ...
{'How many contrasts of interest?';'ncontrast'} , 1);
if strcmpi(button, 'cancel') || isempty(button), return; end
vec = repmat('0 ', 1, params.nconds);
all = [];
for c = 1:condata.ncontrast
tmp = [{{sprintf('Vector for Contrast %d', c); sprintf('con%d', c)}},
{vec},
{{sprintf('Weight for Contrast %d', c); sprintf('con%dw', c)}},
{1}];
all = [all; tmp];
end
[data2, button] = settingsdlg(...
'title', 'Contrast Specification', ...
all{:});
if strcmpi(button, 'cancel') || isempty(button), return; end
con = struct2cell(data2);
params.L = [];
convec = con(1:2:end);
conweight = con(2:2:end);
params.L = [];
params.conWeights = [];
for c = 1:length(conweight)
params.L = [params.L; str2num(convec{c})];
params.conWeights(c) = conweight{c};
end
% =========================================================================
% Save Design Parameters
% =========================================================================
[d, t] = get_timestamp;
outfile = sprintf('design_params_%s_%s.mat', d, t);
fprintf('\nParameters saved in %s\n', fullfile(pwd, outfile));
save(fullfile(pwd, outfile), 'params');
end
function params = defineX(params)
% Pulse Sequence Parameters
params.nslices = 16; % number of time bins in each TR (SPM default = 16)
% Derive Some Additional Parameters %
params.ntrials=sum(params.trialsPerCond); % computes total number of trials
params.scan_length=ceil((params.restBegin + params.restEnd + params.ntrials*(params.meanISI+params.trialDur))/params.TR); % computes total scan length (in TRs)
params.TReff=params.TR/params.nslices; % computes effective TR
%% Verify that valid order can be made with specified max reps %%
fprintf('\nVerifying efficiency of creating valid trial orders with current parameters... ');
tic; % start the timer
order = [];
for i = 1:params.nconds, order = [order; repmat(i, params.trialsPerCond(i), 1)]; end
move_on = 0;
maxtesttime = 5;
while ~move_on
tmp = order(randperm(params.ntrials));
nchunk = getchunks(tmp);
if ~any(nchunk>params.maxRep)
move_on = 1;
elseif toc>=maxtesttime
fprintf('FAIL\n\n');
error('Current parameters prevent efficient creation of valid trial orders. Try increasing Maximum Block Size parameter.');
end
end
order = tmp;
% if toc>=maxtime*60, break, end
fprintf('SUCCESS\n');
% Get ISI distribution %
fprintf('\nCreating base distirbution of ISIs... ');
minISI = params.minISI;
maxISI = params.maxISI;
meanISI = params.meanISI;
TReff = params.TReff;
mu = TReff:TReff:meanISI;
jitSample = zeros(1000,length(mu));
% try
% for s = 1:length(mu), jitSample(:,s) = random(params.distISI, mu(s), 1000, 1); end
% catch
switch lower(params.distISI)
% Thanks to Jeremy Purcell for this fix for older versions of MATLAB
case 'rayleigh'
for s = 1:length(mu), jitSample(:,s) = raylrnd(mu(s), 1000, 1); end
case 'chisquare'
for s = 1:length(mu), jitSample(:,s) = chi2rnd(mu(s), 1000, 1); end
case 'exponential'
for s = 1:length(mu), jitSample(:,s) = exprnd(mu(s), 1000, 1); end
end
% end
jitSample(jitSample<minISI) = NaN;
jitSample(jitSample>maxISI) = NaN;
jitdist = abs(meanISI - nanmean(jitSample));
minIDX = find(jitdist==min(jitdist));
params.jitSample = jitSample(:,minIDX(1));
params.jitSample(isnan(params.jitSample)) = [];
fprintf('SUCCESS\n');
end
function optimizeX(params)
keep = params.keep;
L = params.L;
conWeights = params.conWeights;
gensize = params.gensize;
gensize = gensize - mod(gensize,2); % ensures gensize is divisible by 2
ngen = params.ngen;
maxtime = params.maxtime;
nalpha = max([params.keep ceil(gensize*.01)]);
halfgen = gensize/2;
L(:,end+1) = 0;
L = L';
ncontrasts = size(L,2);
ngenbin = 10;
%% Check Settings %%
if size(L,1)~=params.nconds+1, error('# of columns in contrast matrix ''L'' does not equal # of conditions defined in params'); end
if length(conWeights)~=size(L,2), error('# of contrast weights does not equal # of contrasts'); end
%% Begin Optimization %%
[d, t] = get_timestamp;
fprintf('\nDesign Optimization Started %s on %s', t, d);
fprintf('\n\n\tDESIGN PARAMETERS\n');
disp(params)
%% Start the Progress Monitor %%
progstr = sprintf('Generation %%0%dd/%d', length(num2str(ngen)), ngen);
tic; % start the timer
h = waitbar(1/ngen,sprintf(progstr, 1),'Name','optimizeXGUI Progress', ...
'CreateCancelBtn','setappdata(gcbf,''canceling'',1)');
setappdata(h,'canceling',0)
%% Create First Generation %%
fprintf(sprintf(['\n' progstr ' '], 1));
efficiency = zeros(gensize,1);
order = cell(gensize,1);
jitter = cell(gensize,1);
genbins = floor(gensize/ngenbin:gensize/ngenbin:gensize);
for i = 1:gensize
d = makeX(params);
X=d.X;
X(:,end+1) = 1;
for c = 1:ncontrasts
eff(c) = 1/trace(L(:,c)'*pinv(X'*X)*L(:,c));
end
efficiency(i) = eff*conWeights';
order{i} = d.combined(:,2);
jitter{i} = d.combined(:,5);
if ismember(i,genbins), fprintf('.'), end
end
fprintf(' Max Efficiency = %2.15f', max(efficiency));
maxgeneff(1) = max(efficiency);
%% Loop Over Remaining Generations %%
for g = 2:ngen
fprintf(sprintf(['\n' progstr ' '], g))
waitbar(g/ngen, h, sprintf(progstr, g));
%% Check for Cancel button press
if getappdata(h,'canceling'), delete(h); break; end
%% Grab the Alphas %%
tmp = sortrows([(1:length(efficiency))' efficiency], -2);
fitidx = tmp(1:nalpha,1);
fit.efficiency = efficiency(fitidx);
fit.order = order(fitidx);
fit.jitter = jitter(fitidx);
%% Use the Alphas to Breed %%
cross.efficiency = zeros(halfgen,1);
cross.order = cell(halfgen,1);
cross.jitter = cell(halfgen,1);
genbins = floor(halfgen/(ngenbin/2):halfgen/(ngenbin/2):halfgen);
for i = 1:halfgen
%% Combine Orders %%
conidx = randperm(params.nconds);
orderidx = randperm(length(fit.order));
fixcon = conidx(1);
varcon = conidx(2:end);
calpha = fit.order{orderidx(1)};
mate = fit.order{orderidx(2)};
calpha(ismember(calpha,varcon)) = mate(ismember(mate,varcon));
d=makeX(params, calpha);
X=d.X;
X(:,end+1) = 1;
for c = 1:ncontrasts
eff(c) = 1/trace(L(:,c)'*pinv(X'*X)*L(:,c));
end
cross.efficiency(i) = eff*conWeights';
cross.order{i} = d.combined(:,2);
cross.jitter{i} = d.combined(:,5);
if ismember(i,genbins), fprintf('.'), end
end
%% Check for Cancel button press
if getappdata(h,'canceling'), delete(h); break; end
%% Introduce Some Nasty Mutants %%
if g>2 && maxgeneff(g-1)==maxgeneff(g-2)
mutsize = gensize;
else
mutsize = halfgen;
end
mut.efficiency = zeros(mutsize,1);
mut.order = cell(mutsize,1);
mut.jitter = cell(mutsize,1);
genbins = floor(mutsize/(ngenbin/2):mutsize/(ngenbin/2):mutsize);
for i = 1:mutsize
d=makeX(params);
X=d.X;
X(:,end+1) = 1;
for c = 1:ncontrasts
eff(c) = 1/trace(L(:,c)'*pinv(X'*X)*L(:,c));
end
mut.efficiency(i) = eff*conWeights';
mut.order{i} = d.combined(:,2);
mut.jitter{i} = d.combined(:,5);
if ismember(i,genbins), fprintf('.'), end
end
%% Combine this Genertation and Compute Max Efficiency %%
efficiency = [fit.efficiency; cross.efficiency; mut.efficiency];
order = [fit.order; cross.order; mut.order];
jitter = [fit.jitter; cross.jitter; mut.jitter];
fprintf(' Max Efficiency = %2.15f', max(efficiency));
maxgeneff(g) = max(efficiency);
%% Break if Over Time %%
if toc>=maxtime*60, break, end
%% Check for Cancel button press
if getappdata(h,'canceling'), delete(h); break; end
end
delete(h);
%% Save Best Designs %%
[d, t] = get_timestamp;
outdir = sprintf('best_designs_%s_%s', d, t); mkdir(outdir);
fprintf('\n\nSaving %d best designs to: %s\n', keep, fullfile(pwd, outdir));
tmp = sortrows([(1:length(efficiency))' efficiency], -2);
fitidx = tmp(1:keep,1);
best.efficiency = efficiency(fitidx);
best.order = order(fitidx);
best.jitter = jitter(fitidx);
design = cell(keep,1);
for i = 1:keep
design{i}=breedX(params, best.order{i}, best.jitter{i});
fname = [outdir filesep 'design' num2str(i) '.txt'];
dlmwrite(fname, design{i}.combined, 'delimiter', '\t')
fname2 = [outdir filesep 'design' num2str(i) '.csv'];
cc = [{'Trial' 'Condition' 'Onset' 'Duration' 'ISI'}; num2cell(design{i}.combined)];
writedesign(cc, fname2);
end
save([outdir filesep 'designinfo.mat'], 'design', 'params');
fprintf('Finished in %2.2f minutes at %s on %s\n\n', toc/60, t, d);
%% Visualize the Best Design
hfig = figure('color', 'white', 'units', 'norm', 'pos', [.3 .3 .3 .4]);
winner = scalemat(design{1}.X);
winner = [winner ones(size(winner,1), 1)];
h = imagesc(winner, 'parent', gca); colormap('gray');
set(gca, 'FontName', 'Arial', 'FontSize', 18, 'XTick', 1:size(winner,2));
xticklabel = strcat({'cond'}, get(gca, 'XTickLabel'));
xticklabel{end} = 'constant';
set(gca, 'xticklabel', xticklabel);
ylabel('TR', 'fontsize', 18, 'fontweight', 'bold');
title('The "Best" Design Matrix', 'fontsize', ceil(18*1.10), 'fontweight', 'bold');
end
function design = makeX(params, order, verbose)
if nargin < 3, verbose = 0; end
if nargin < 2
makeorder = 1;
elseif isempty(order)
makeorder = 1;
else
makeorder = 0;
end
%-----------------------------------------------------------------
% Get a pseudoexponential distribution of ISIs
%-----------------------------------------------------------------
goodJit=0;
while goodJit==0
jitters=randsample(params.jitSample,params.ntrials-1,1);
if mean(jitters) < params.meanISI+params.TReff && mean(jitters) > params.meanISI-params.TReff
goodJit=1;
end
end
%-----------------------------------------------------------------
% Determine stimulus onset times
%-----------------------------------------------------------------
onset=zeros(1,params.ntrials);
onset(1)=params.restBegin;
for t=2:params.ntrials,
onset(t)=onset(t-1) + params.trialDur + jitters(t-1);
end;
jitters(end+1)=params.restEnd;
%-----------------------------------------------------------------
% Make some trial orders
%-----------------------------------------------------------------
if makeorder
if params.counterBalanceInterval
order = make_order(params.trialsPerCond, params.maxRep, params.counterBalanceInterval);
else
order = [];
for i = 1:params.nconds, order = [order; repmat(i, params.trialsPerCond(i), 1)]; end
move_on = 0;
while ~move_on
tmp = order(randperm(params.ntrials));
nchunk = getchunks(tmp);
if ~any(nchunk>params.maxRep), move_on = 1; end
end
order = tmp;
end
end
%------------------------------------------------------------------------
% Create the design matrix (oversample the HRF depending on effective TR)
%------------------------------------------------------------------------
cond=order;
oversamp_rate=params.TR/params.TReff;
dmlength=params.scan_length*oversamp_rate;
oversamp_onset=(onset/params.TR)*oversamp_rate;
hrf=spm_hrf(params.TReff);
desmtx=zeros(dmlength,params.nconds);
for c=1:params.nconds
r=zeros(1,dmlength);
cond_trials= cond==c;
cond_ons=fix(oversamp_onset(cond_trials))+1;
r(cond_ons)=1;
cr=conv(r,hrf);
desmtx(:,c)=cr(1:dmlength)';
onsets{c}=onset(cond==c); % onsets in actual TR timescale
end;
% sample the design matrix back into TR timescale
desmtx=desmtx(1:oversamp_rate:dmlength,:);
%------------------------------------------------------------------------
% Filter the design matrix
%------------------------------------------------------------------------
K.RT = params.TR;
K.HParam = params.hpf;
K.row = 1:length(desmtx);
K = spm_filter(K);
for c=1:params.nconds
desmtx(:,c)=spm_filter(K,desmtx(:,c));
end
%------------------------------------------------------------------------
% Save the design matrix
%------------------------------------------------------------------------
design.X = desmtx;
design.combined=zeros(params.ntrials,5);
design.combined(:,1)=1:params.ntrials;
design.combined(:,2)=cond;
design.combined(:,3)=onset;
design.combined(:,4)=repmat(params.trialDur,params.ntrials,1);
design.combined(:,5)=jitters;
design.duration=(params.scan_length*params.TR)/60;
end
function design = breedX(params, order, jitters)
% USAGE: design = breedX(params, order, jitters)
if nargin<3, display('USAGE: design = breedX(params, order, jitters)'), return, end
%-----------------------------------------------------------------
% Determine stimulus onset times
%-----------------------------------------------------------------
onset=zeros(1,params.ntrials);
onset(1)=params.restBegin;
for t=2:params.ntrials,
onset(t)=onset(t-1) + params.trialDur + jitters(t-1);
end;
%------------------------------------------------------------------------
% Create the design matrix (oversample the HRF depending on effective TR)
%------------------------------------------------------------------------
cond=order;
oversamp_rate=params.TR/params.TReff;
dmlength=params.scan_length*oversamp_rate;
oversamp_onset=(onset/params.TR)*oversamp_rate;
hrf=spm_hrf(params.TReff);
desmtx=zeros(dmlength,params.nconds);
for c=1:params.nconds
r=zeros(1,dmlength);
cond_trials= cond==c;
cond_ons=fix(oversamp_onset(cond_trials))+1;
r(cond_ons)=1;
cr=conv(r,hrf);
desmtx(:,c)=cr(1:dmlength)';
onsets{c}=onset(cond==c); % onsets in actual TR timescale
end;
% sample the design matrix back into TR timescale
desmtx=desmtx(1:oversamp_rate:dmlength,:);
%------------------------------------------------------------------------
% Filter the design matrix
%------------------------------------------------------------------------
K.RT = params.TR;
K.HParam = params.hpf;
K.row = 1:length(desmtx);
K = spm_filter(K);
for c=1:params.nconds
desmtx(:,c)=spm_filter(K,desmtx(:,c));
end
%------------------------------------------------------------------------
% Save the design matrix
%------------------------------------------------------------------------
design.X = desmtx;
design.combined=zeros(params.ntrials,5);
design.combined(:,1)=1:params.ntrials;
design.combined(:,2)=cond;
design.combined(:,3)=onset;
design.combined(:,4)=repmat(params.trialDur,params.ntrials,1);
design.combined(:,5)=jitters;
design.duration=(params.scan_length*params.TR)/60;
end
% =========================================================================
% *
% * SUBFUNCTIONS (ORDERING)
% *
% =========================================================================
function ordervec = make_order(condtrials, maxrep, cbinterval, cborder)
% MAKE_ORDER Make trial order w/optional counterbalancing
%
% USAGE: [ordervec, orderbinmat] = make_order(condtrials, maxrep, [cbinterval], [cborder])
%
% ARGUMENTS
% condtrials = vector whose length corresponds to the number of
% conditions and whose values correspond to the number of trials per
% condition. For example, [10 10 15] indicates that there are 3
% conditions, with conditions 1 and 2 featuring 10 trials and
% condition 3 featuring 15 trials.
% maxrep = scalar that indicate the maximum number of times a trial
% is allowed to be repeated on consecutive trials. Use to control the
% extent to which trials from the same condition can be blocked.
% cbinterval = counterbalancing across equal intervals
% 0 or 1 (default) = none
% > 1 = will divide into equal intervals of the specified size, e.g.,
% 2 will counterbalance across first and second half
% 4 will counterbalance across quartiles
% cborder = order of m-seq counterbalancing (0 for none)
%
% --------------------------------------- Copyright (C) 2014 ---------------------------------------
% Author: Bob Spunt
% Affilitation: Caltech
% Email: spunt@caltech.edu
%
% $Revision Date: Aug_20_2014
if nargin < 2, error('USAGE: [ordervec, orderbinmat] = make_order(condtrials, maxrep, [cbinterval], [cborder])'); end
if nargin < 3, cbinterval = 0; end
if nargin < 4, cborder = 0; end
if maxrep==1, error('maxrep must be greater than 1'); end
ncond = length(condtrials);
ntrials = sum(condtrials);
ntrialsmax = ncond*max(condtrials);
if cbinterval > 1
qt = 0:(1/cbinterval):1;
goodn = [floor(condtrials/cbinterval); ceil(condtrials/cbinterval)];
qidx1 = ceil(quantile(1:ntrials, qt(1:end-1)));
qidx2 = [qidx1(2:end)-1 ntrials];
end
if cborder
seqrep = 0;
tmp = 0;
while length(tmp) < ntrialsmax
seqrep = seqrep + 1;
tmp = carryoverCounterbalance(ncond, cborder, seqrep);
end
move_on = 0;
while ~move_on
move_on = 1;
seq = carryoverCounterbalance(ncond, cborder, seqrep)';
for i = 1:ncond
condseq = find(seq==i);
seq(condseq(condtrials(i)+1:end)) = NaN;
end
tmp = seq(~isnan(seq));
if any(getchunks(tmp) > maxrep), move_on = 0; continue; end
if cbinterval > 1
orderbinmat = repmat(tmp, 1, ncond)==repmat(1:ncond, ntrials, 1);
for i = 1:cbinterval
if ~all(sum(repmat(sum(orderbinmat(qidx1(i):qidx2(i), :)), size(goodn, 1), 1)==goodn))
move_on = 0; continue; end
end
end
end
ordervec = tmp;
else
ordervec = [];
for i = 1:ncond, ordervec = [ordervec; repmat(i, condtrials(i), 1)]; end
move_on = 0;
while ~move_on
move_on = 1;
tmp = ordervec(randperm(ntrials));
if any(getchunks(tmp) > maxrep), move_on = 0; continue; end
if cbinterval > 1
orderbinmat = repmat(tmp, 1, ncond)==repmat(1:ncond, ntrials, 1);
for i = 1:cbinterval
if ~all(sum(repmat(sum(orderbinmat(qidx1(i):qidx2(i), :)), size(goodn, 1), 1)==goodn))
move_on = 0; continue; end
end
end
end
ordervec = tmp;
end
end
function order = make_order_old(condtrials, maxrep, cborder)
if nargin < 3, error('USAGE: order = make_order(condtrials, maxrep, cborder)'); end
if maxrep==1, error('maxrep must be greater than 1'); end
ncond = length(condtrials);
ntrials = sum(condtrials);
seqrep = ceil(ntrials/length(carryoverCounterbalance(ncond, cborder, 1)));
move_on = 0;
while ~move_on
seq = carryoverCounterbalance(ncond, cborder, seqrep);
for i = 1:ncond
condseq = find(seq==i);
seq(condseq(condtrials(i)+1:end)) = NaN;
end
tmp = seq(~isnan(seq));
if ~any(getchunks(tmp) > maxrep), move_on = 1; end
end
order = tmp;
end
function condSequence = carryoverCounterbalance(numConds,cbOrder,reps)
% CARRYOVERCOUNTERBALANCE
%
% USAGE: condSequence = carryoverCounterbalance(numConds,cbOrder,reps)
%
% This is a minimally modified version of software written by Joseph Brooks
% and described in the following paper which should be cited if use
% counterbalancing:
%
% Brooks, J.L. (2012). Counterbalancing for serial order carryover effects
% in experimental condition orders. Psychological Methods. Vol 17(4), Dec
% 2012, 600-614. doi.org/10.1037/a0029310
%
%__________________________________________________________________________
% Creates a serial-order counterbalanced sequence of conditions, i.e. it
% ensures that every condition is preceded equally often by every other
% condition (first-order) or every nTuple of length n for nth-order
% counterbalancing. This can be done for any number of conditions. Uses
% Kandel et al. 1996 method to determine Euler path (see accessory
% functions included below).
%
% OUTPUT: -condSequence: the counterbalanced sequence. Length will depend
% on parameters. Output as an array of numbers
%
% INPUTS:
% -numConds: N unique conditions. Must be a positive integer
%
% -cbOrder: depth/order of counterbalancing. Must be a positive integer.
% First-order counterbalancing ensures that every condition is preceded
% by every other condition equally often. nth-order counterbalancing
% ensures that every condition is preceded equally often by every
% n-length sequence (nTuple) of conditions. e.g. 2nd-order
% counterbalancing would ensure that condition 1 is preceded equally
% often by each set of 2 conditions, 1 1, 1 2, 2 1, 2 2.
%
% -reps: The number of times each condition is preceded by every other
% condition (first-order) or every nTuple (nth-order). Must be a positive
% integer. This can be used to increase the total number of trials while
% maintaining counterbalancing.
%
% -omitSelfAdjacencies: 0 = include condition self-adjacencies (e.g. 2 2
% 1 2 1 1 2); 1 = exclude condition self-adjacencies (e.g. 1 2 1 2 1 2)
%
% VERSION: 1.0.01.03.2012
% Joseph Brooks, UCL Institute of Cognitive Neuroscience brooks.jl@gmail.com
%
omitSelfAdjacencies = 0;
if nargin<1, error('USAGE: condSequence = carryoverCounterbalance(numConds,cbOrder,rep)'); end
if nargin<2, cbOrder = 1; end
if nargin<3, reps = 2; end
if reps < 1, error('ERROR: reps parameter must be > 0'); end
if cbOrder < 1, error('ERROR: cbOrder parameter must be > 0'); end
if numConds < 1, error('ERROR: numConds parameter must be > 0'); end
if omitSelfAdjacencies < 0 || omitSelfAdjacencies > 1,
error('ERROR: omitSelfAdjacencies parameter must be 0 or 1');
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CONSTRUCT ADJACENCY MATRIX FOR GRAPH WITH APPROPRIATE TEMPORAL
% RELATIONSHIPS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if cbOrder == 1
%%%Construct adjacency matrix and condition list for cbOrder == 1
%%%Fully-connected graph
adjacencyMatrix = ones(numConds,numConds)*reps;
nTuples = [1:numConds];
if omitSelfAdjacencies
adjacencyMatrix = adjacencyMatrix - adjacencyMatrix.*eye(size(adjacencyMatrix));
end;
else
%%%CONSTRUCT N-TUPLE CORRESPONDING TO EACH NODE WHEN cbOrder > 1
nTuples = nTuples_brooks(numConds,cbOrder);
%If indicated in input parameters, remove nodes with self-adjacencies by putting zeros in
%corresponding columns and rows of adjacency matrix
if omitSelfAdjacencies
nTuples = nTuples_removeSelfAdjacencies_brooks(nTuples);
end
%%%Construct adjacency matrix by connecting only those nTuples
%%%which share the (n-1)-length-beginning/(n-1)-length-end of their sequences
adjacencyMatrix = zeros(size(nTuples,1),size(nTuples,1));
for tminus1 = 1:size(adjacencyMatrix,1)
for t = 1:size(adjacencyMatrix,2)
if nTuples(tminus1,2:size(nTuples,2)) == nTuples(t,1:size(nTuples,2)-1)
adjacencyMatrix(tminus1,t) = 1;
end
end
end
%%%Duplicate edges based on number of reps requested
adjacencyMatrix = adjacencyMatrix*reps;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%FIND EULER PATH THROUGH GRAPH SPECIFIED BY ADJACENCY MATRIX%%%Uses Kandel et al 1996 method for Euler Circuit
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nodeSequence = eulerPathKandelMethod_carryoverCounterbalance(adjacencyMatrix);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%CONSTRUCT CONDITION SEQUENCE FROM NODE SEQUENCE.
%%%-For first-order counterbalancing node sequence = cond sequence
%%%-For higher-order counterbalancing, overlapping parts of sequence
%%%must be considered and only one additional condition from the
%%%nTuple will be added for all nodes beyond the first
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if cbOrder == 1
condSequence = nodeSequence;
else
condSequence = nTuples(nodeSequence(1),:);
for i = 2:length(nodeSequence)
condSequence = [condSequence nTuples(nodeSequence(i),size(nTuples,2))];
end
end;
end
function [seq] = eulerPathKandelMethod_carryoverCounterbalance(adjacencyMatrix)
% Finds an Euler circuit through a graph, G, and returns sequence of nodes
% associated with that path. Based on method by Kandel, et al. (1996).
% The circuit/path created is uniformly drawn from the set of all possible
% Eulerian circuits of G.
%
% USAGE:
% [seq] = eulerPathKandelMethod(adjacencyMatrix)
% OUTPUT: seq is a ordered sequence of nodes corresponding to the circuit
% through G.
%
% INPUT: adjacency matrix
% Option1: If entering a single integer, n, then program assumes a
% fully-connected digraph with n nodes and 1 instance of the edge
% between each node.
% Option2: If entering a matrix then the matrix will be
% treated as the adjacency matrix of the graph, G. Must be an nxn matrix for
% a graph with n nodes Dimension1 = represents the
% node sending the directed edge, e.g. G(1,2) is an entry
% representing a directed edge from node 1 to node 2.
%
% - adjacency matrix must be consistent with an Eulerian graph,
% i.e., all nodes have in-degree = out-degree.
% Semi-Eulerian graphs are not valid for this program.
% - all entries in the adjacency matrix must be >= 0. Values
% larger than one indicate duplicated instances of an edge between
% the nodes.
%
% References:
% Kandel, D., Matias, Y., Unger, R., & Winkler, P. (1996). Shuffling
% biological sequences. Discrete Applied Mathematics, 71(1-3), 171-185:
%
% VERSION: 1.0.01.03.2012
% by Joseph Brooks, UCL Institute of Cognitive Neuroscience
% brooks.jl@gmail.com%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% CHECK INPUT & GIVE ERROR MESSAGES IF NOT APPROPRIATE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%If input is integer then create adjacency matrix for fully connected
%%%digraph with one instance of each edge
if size(adjacencyMatrix,1) == 1
adjacencyMatrix = ones(adjacencyMatrix,adjacencyMatrix);
end
%%%Is matrix square?
if size(adjacencyMatrix,1) ~= size(adjacencyMatrix,2)
error('ERROR: Adjacency matrix is not square');
end;
%%%Is matrix 2D?
if length(size(adjacencyMatrix)) > 2
error('ERROR: Adjacency matrix should have only 2 dimensions');
end;
%%% Is graph Eulerian? Error if not. Error if semi-Eulerian
for i = 1:size(adjacencyMatrix,1)
if sum(adjacencyMatrix(:,i)) ~= sum(adjacencyMatrix(i,:))
error('ERROR: non_Eulerian graph specfied. In-degree must equal out-degree at every vertex');
end;
end
%%% Does each vertex have at least one edge?
for i = 1:size(adjacencyMatrix,1)
if sum(adjacencyMatrix(:,i)) == 0 && sum(adjacencyMatrix(i,:)) == 0
error('ERROR: Disconnected graph...at least one vertex has no edges.');
end;
end
%%%Are all values >= 0
if min(min(adjacencyMatrix)) < 0
error('ERROR: Adjacency matrix contains value < 0');
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% START COMPUTATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Uniformly draw arborescence from the set of all possible for the input graph
%%%arb output is adjacency matrix for an arborescence of the input
arb = arborescenceKandel_carryoverCounterbalance(adjacencyMatrix);
%The following matrix indicates which out edges are NOT included in the
%arborescence by subtracting the arborescence from the adjacency matrix.
%These edges need to be randomly permuted in the outEdges list
remainingEdges = adjacencyMatrix - arb;
%%%For each node create a list of outgoing edges and randomly order them
%%%except for any edge that is part of the arborescence. Any outgoing edges
%%%that are in the arborescence should be last in the order
outEdgesOrder = cell(size(adjacencyMatrix,1),1);for o = 1:size(adjacencyMatrix,1)
%%%In the cell corresponding to node n, list all out edges from this
%%%node including duplicate out edges as muliple instances in the list.
for i = 1:size(adjacencyMatrix,2)
for r = 1:remainingEdges(o,i)
outEdgesOrder{o} = [outEdgesOrder{o} i];
end
end
%%%Randomly shuffle the out edges for this node
outEdgesOrder{o} = outEdgesOrder{o}(randperm(length(outEdgesOrder{o})));
%%%Add arborescence edge to end of list if it exists
if sum(arb(o,:)) > 0
outEdgesOrder{o} = [outEdgesOrder{o} find(arb(o,:) == 1)];
end
end
%%%Set first node in sequence to the root of the arboresence (node
%%%toward which all edges in arb are pointing)
seq = find(sum(arb,2) == 0);
%%% Generate sequence of nodes by starting at the root node of the arb
%%% matrix and take the out edge from that node with the lowest number.
%%% Add the target node of this out edge to the sequence and Remove the out
%%% edge from the list for the source node. Repeat treating the target node
%%% as the new source node until all edges have been traversed.
while sum(cellfun('length',outEdgesOrder)) > 0
seq = [seq outEdgesOrder{seq(end)}(1)];
if length(outEdgesOrder{seq(end-1)}) > 1
outEdgesOrder{seq(end-1)} = outEdgesOrder{seq(end-1)}(2:end);
else
outEdgesOrder{seq(end-1)} = [];
end
end
end
function arb = arborescenceKandel_carryoverCounterbalance(adjacencyMatrix)
%%% FIND UNIFORMLY-DRAWN ARBORESCENCE for the graph specified by the adjacency matrix
%%% in the input. Do this by performing a backward random walk
%%% on the graph...based on Kandel et al. (1996)
%%%
%%% INPUT: Adjacency Matrix of the graph. A square matrix of values >= 0
%%% where a positive integer in Aij indicates an edge from vertex i to
%%% vertex j. The size of each dimension of the matrix is equal to the
%%% number of vertices of the graph
%%%
%%% OUTPUT: Adjacency Matrix of the arborescence.
%%%
%%% KANDEL ALGORITHM
%%% (1) start at random node
%%% (2) randomly choose an edge to cross (loops which connect back to the same node can be ignored)
%%% (3) if this edge leads to a node that has not been visited then
%%% this edge is part of the arborescence. For the edge's origin node,
%%% it should be marked as the last edge to be traversed out of the
%%% node when creating the sequence.
%%% in outEdgeOrder matrix to make it the last traversed, i.e. size-of-matrix+1%%% (4) select an edge out and move to another node.
%%% (5) repeated 2-4 until each node has been visited at least once.
%%%
%%% VERSION: 1.0.01.03.2012
%%% by Joseph Brooks, UCL Institute of Cognitive Neuroscience
%%% brooks.jl@gmail.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%CHECK INPUT%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Is matrix square?
if size(adjacencyMatrix,1) ~= size(adjacencyMatrix,2)
error('ERROR: Adjacency matrix is not square');
end;
%%%Is matrix 2D?
if length(size(adjacencyMatrix)) > 2
error('ERROR: Adjacency matrix should have only 2 dimensions');
end;
%%% Does each vertex have at least one edge?
for i = 1:size(adjacencyMatrix,1)
if sum(adjacencyMatrix(:,i)) == 0 && sum(adjacencyMatrix(i,:)) == 0
error('ERROR: At least one vertex is disconnected (i.e. has no edges)');
end;
end
%%%Are all values >= 0
if min(min(adjacencyMatrix)) < 0
error('ERROR: Adjacency matrix contains value < 0');
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%COMPUTE ARBORESCENCE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Create empty adjacency matrix to represent output arborescence
arb = zeros(size(adjacencyMatrix));
%%%Simplify input adjacency matrix by removing duplicate edges. The result
%%%will treat duplicate edges as equivalent.
adjacencyMatrix = sign(adjacencyMatrix);
%%%Choose a random starting vertex and add it to the list of nodes visited
currentNode = randi(size(adjacencyMatrix,1));
nodesVisited = [currentNode];
while length(unique(nodesVisited)) < size(adjacencyMatrix,1)
%%%Find all edges INTO the current node...edges into the node are
%%%considered because this is a BACKWARD walk
availableSourceNodes = find(adjacencyMatrix(:,currentNode) > 0);
%%%Randomly choose one of the edges into of the node and designate as
%%%source
selectedSource = availableSourceNodes(randi(length(availableSourceNodes)));
%%%If this is the first visit to the source node, then mark the edge as
%%%part of the arborescence
if sum([nodesVisited == selectedSource]) == 0
arb(selectedSource,currentNode) = 1; end
%%%Add target node to list of visited nodes
nodesVisited = [nodesVisited,selectedSource];
currentNode = selectedSource;
end
end
function result = nTuples_brooks(numItems,n)
%
% Create all possible nTuples of length n from a list of items of length =
% numItems
%
% OUTPUT:
% -result: matrix where each row corresponds to an nTuple of length n. Size
% of matrix will be numItems^n x n
%
% INPUT:
% -numItems: this is the number of items that are possible in each position
% of each nTuple.
% -n: this is the length of each nTuple.
%
% EXAMPLE: For all of the nTuples of length 2 of 3 items, use nTuples(3,2).
% The result of this computation is:
% 1 1
% 1 2
% 1 3
% 2 1
% 2 2
% 2 3
% 3 1
% 3 2
% 3 3
%
% VERSION: 1.0.05.03.2012
% by Joseph Brooks, UCL Institute of Cognitive Neuroscience
% brooks.jl@gmail.com
result = zeros(numItems^n,n);
for v = 1:numItems^n
for i = 1:n
result(v,i) = mod(ceil(v/numItems^(i-1)),numItems)+1;