Click here to download this file

Click here to download a scripted example

function [Out] = SamsOptseq(In)
% "Sam's OptSeq"
% By Sam Berens (S.Berens@sussex.ac.uk).
%
% [Out] = SamsOptSeq(In)
%
% INPUT: "In" should be a 1x1 structure with the following fields:
%   - "In.nEventReps" : a nx1 numeric denoting the number of repetitions...
%      for each of n experimental event types;
%   - "In.ContrastMats" : an cx1 cell array of axn matrices. This specifies...
%     "c" contrasts where "a" is the number of rows in any one contrast,...
%     and "n" is the number of experimental conditions (or regressors), as...
%     specified in "In.nEventReps";
%   - "In.ContrastNames" : an cx1 cell array of strings specifying the names...
%     of each contrast in "In.ContrastMat";
%   - "In.nNulls" : a 1x1 numeric denoting the number of null events...
%     (unmodelled) to be induced;
%   - "In.ITI" : a 1x1 numeric denoting the inter-trial-interval (ITI) in...
%     seconds for both experimental and null events;
%   - "In.nScans" : a 1x1 numeric denoting the number of null scans to be...
%     acquired;
%   - "In.TR" : a 1x1 numeric denoting the repetition time (TR) of the...
%     scanning sequence in seconds;
%   - "In.nIterations" : a 1x1 numeric denoting the number of iterations...
%     that the algorithm should search through;
%   - "In.CutoffFq_Sc" : [Optional] a 1x1 numeric denoting the high-pass...
%     cuff-off frequency (in seconds) used to filter the fMRI dataset...
%     (default = 128s);
%
% OUTPUT: an nx1 structure that lists each of the tested stimulus...
% order permutations alongside the efficiency statistics for each of...
% the specified contrasts. Efficiency statistics are monotonically (but
% perhaps not parametrically) related to statistical power. They are
% relative measures that depend on the scaling of the design matrix and
% the scaling of the contrasts that were specified.
%   For more information, see...
%       http://imaging.mrc-cbu.cam.ac.uk/imaging/DesignEfficiency.
%   Thanks to Rik Henson for his insightful fMRI pages!

% --- Error checking:
StimulationDuration = (sum(In.nEventReps,1) + In.nNulls) * In.ITI;
ScanDuration = In.nScans * In.TR;
if StimulationDuration > ScanDuration
    warning('The duration of experimental stimulation (%f seconds)%cis greater than the length of the scanning sequence (%f seconds).',...
        StimulationDuration,10,ScanDuration);
end
for i_Contrasts = 1:1:size(In.ContrastMats,1)
    if size(In.ContrastMats{i_Contrasts,1},2) ~= size(In.nEventReps,1)
        error('Regressor count mismatch: size(In.nEventReps,1) ~= size(In.ContrastMats{%i,1},2).',i_Contrasts);
    end
end
if size(In.ContrastMats,1) ~= size(In.ContrastNames,1)
    error('Contrast count mismatch: size(In.ContrastMats,1) ~= size(In.ContrastNames,1).');
end

global CutoffFq_Sc;
if ~isfield(In,'CutoffFq_Sc')
    CutoffFq_Sc = 128;
else
    CutoffFq_Sc = In.CutoffFq_Sc;
end

% --- Initialise values:
eHrf_MS = MakeEffectiveHrf(0.001);
eHrf_MS = eHrf_MS ./ max(eHrf_MS,[],1);
nAll = zeros((size(In.nEventReps,1)+1),1);
nAll(1,1) = In.nNulls;
nAll(2:end,1) = In.nEventReps;
Out = struct;
dispstat('','init');

% --- Iteration loop:
for i_Iteration = 1:1:In.nIterations

    % Make StimOrder:
    StimOrder = nan(sum(nAll,1),2);
    i_StartWrite = 1;
    for i_TrialType = 1:1:size(nAll,1)
        i_EndWrite = (i_StartWrite + nAll(i_TrialType,1)) - 1;
        StimOrder(i_StartWrite:1:i_EndWrite,2) = i_TrialType - 1;
        i_StartWrite = i_EndWrite + 1;
    end
    StimOrder(:,1) = rand(size(StimOrder,1),1);
    StimOrder = sortrows(StimOrder,1);
    StimOrder = StimOrder(:,2);
    Out(i_Iteration,1).StimOrder = StimOrder;

    % Make DesignMatrix_SC_PostConv:
    DesignMatrix_MS_PreConv = zeros((In.nScans*1000*In.TR),size(In.nEventReps,1));
    ScanSamples_MS = mod((0:1:(size(DesignMatrix_MS_PreConv,1)-1))',(1000*In.TR)) == 0;
    i_Event = 0;
    for i_DesignMatrix_MS_PreConv = 1:(In.ITI*1000):(In.ITI*1000*size(StimOrder,1))
        i_Event = i_Event + 1;
        if StimOrder(i_Event,1) ~= 0
            DesignMatrix_MS_PreConv(i_DesignMatrix_MS_PreConv,StimOrder(i_Event,1)) = 1;
        end
    end
    DesignMatrix_MS_PostConv = nan(size(DesignMatrix_MS_PreConv));
    for j_DesignMatrix_MS_PostConv = 1:1:size(DesignMatrix_MS_PostConv,2)
        RegressorTs = conv(DesignMatrix_MS_PreConv(:,j_DesignMatrix_MS_PostConv),eHrf_MS);
        DesignMatrix_MS_PostConv(:,j_DesignMatrix_MS_PostConv) =...
            RegressorTs(1:size(DesignMatrix_MS_PostConv,1),1);
    end
    i_DesignMatrix_SC_PostConv = 0;
    DesignMatrix_SC_PostConv = nan(i_DesignMatrix_SC_PostConv,size(In.nEventReps,1));
    for i_ScanSamples_MS = 1:1:size(ScanSamples_MS,1)
        if ScanSamples_MS(i_ScanSamples_MS,1)
            i_DesignMatrix_SC_PostConv = i_DesignMatrix_SC_PostConv + 1;
            DesignMatrix_SC_PostConv(i_DesignMatrix_SC_PostConv,:) = DesignMatrix_MS_PostConv(i_ScanSamples_MS,:);
        end
    end

    % Compute efficiency for each contrast:
	for i_Contrast = 1:1:size(In.ContrastMats,1)
		ContrastMat = In.ContrastMats{i_Contrast,1};
		Efficiency = trace(ContrastMat*inv(DesignMatrix_SC_PostConv'*DesignMatrix_SC_PostConv)*ContrastMat')^-1; %#ok<MINV,NASGU>
		eval(sprintf('Out(i_Iteration,1).%s = Efficiency;',In.ContrastNames{i_Contrast,1}));
	end
    dispstat(sprintf('Finished iteration %i of %i',i_Iteration,In.nIterations));
end
disp('Done!')
return

function [eHrf_tFunction] = MakeEffectiveHrf(Tr)
global CutoffFq_Sc;
Sr = 1 / Tr; % Sample rate
Hrf_tFunction = spm_hrf(Tr);
Hrf_fFunction = fft(Hrf_tFunction);
BinWidth_Hz = Sr / size(Hrf_fFunction,1);
CutoffFq_Hz = 1 / CutoffFq_Sc;
NumberOfBinsToCut = CutoffFq_Hz / BinWidth_Hz;
eHrf_fFunction = Hrf_fFunction;
eHrf_fFunction(1,1) = Hrf_fFunction(1,1) * (1 - NumberOfBinsToCut);
eHrf_tFunction = ifft(eHrf_fFunction);
return

function dispstat(TXT,varargin)
% Prints overwritable message to the command line. If you dont want to keep
% this message, call dispstat function with option 'keepthis'. If you want to
% keep the previous message, use option 'keepprev'. First argument must be
% the message.
% IMPORTANT! In the firt call, option 'init' must be used for initialization purposes.
% Options:
%     'init'      this must be called in the begining. Otherwise, it can overwrite the previous outputs on the command line.
%     'keepthis'    the message will be persistent, wont be overwritable,
%     'keepprev'  the previous message wont be overwritten. New message will start from next line,
%     'timestamp' current time hh:mm:ss will be appended to the begining of the message.
% Example:
%   clc;
%   fprintf('12345677890\n');
%   dispstat('','init')      %Initialization. Does not print anything.
%   dispstat('Time stamp will be written over this text.'); % First output
%   dispstat('is current time.','timestamp','keepthis'); % Overwrites the previous output but this output wont be overwritten.
%   dispstat(sprintf('*********\nDeveloped by %s\n*********','Kasim')); % does not overwrites the previous output
%   dispstat('','timestamp','keepprev','keepthis'); % does not overwrites the previous output
%   dispstat('this wont be overwriten','keepthis');
%   dispstat('dummy dummy dummy');
%   dispstat('final stat');
% % Output:
%     12345677890
%     15:15:34 is current time.
%     *********
%     Developed by Kasim
%     *********
%     15:15:34
%     this wont be overwriten
%     final stat

% **********
% **** Options
keepthis = 0; % option for not overwriting
keepprev = 0;
timestamp = 0; % time stamp option
init = 0; % is it initialization step?
if ~isstr(TXT) %#ok<DISSTR>
    return
end
persistent prevCharCnt;
if isempty(prevCharCnt)
    prevCharCnt = 0;
end
if nargin == 0
    return
elseif nargin > 1
    for i = 2:nargin
        eval([varargin{i-1} '=1;']);
    end
end
if init == 1
    prevCharCnt = 0;
    return;
end
if isempty(TXT) && timestamp == 0
    return
end
if timestamp == 1
    c = clock; % [year month day hour minute seconds]
    txtTimeStamp = sprintf('%02d:%02d:%02d ',c(4),c(5),round(c(6)));
else
    txtTimeStamp = '';
end
if keepprev == 1
    prevCharCnt = 0;
end
% *************** Make safe for fprintf, replace control charachters
TXT = strrep(TXT,'%','%%');
TXT = strrep(TXT,'\','\\');
% *************** Print
TXT = [txtTimeStamp TXT '\n'];
fprintf([repmat('\b',1, prevCharCnt) TXT]);
nof_extra = length(strfind(TXT,'%%'));
nof_extra = nof_extra + length(strfind(TXT,'\\'));
nof_extra = nof_extra + length(strfind(TXT,'\n'));
prevCharCnt = length(TXT) - nof_extra; %-1 is for \n
if keepthis == 1
    prevCharCnt = 0;
end
return

Sam Berens