From 47290384132005ce60c2cf00bfab115f0d2aa288 Mon Sep 17 00:00:00 2001 From: Daan Wesselink Date: Wed, 6 Sep 2017 13:23:45 +0100 Subject: [PATCH 1/7] What do you think? --- Hello/test.m | 1 + 1 file changed, 1 insertion(+) create mode 100644 Hello/test.m diff --git a/Hello/test.m b/Hello/test.m new file mode 100644 index 0000000..5ab2f8a --- /dev/null +++ b/Hello/test.m @@ -0,0 +1 @@ +Hello \ No newline at end of file From 67ee8a23697510f312a17ac7302bebdfd3b14268 Mon Sep 17 00:00:00 2001 From: Roni Maimon Date: Wed, 6 Sep 2017 14:50:11 +0100 Subject: [PATCH 2/7] setting up the fsl package --- +rsa/+fsl/getDataFromFSL.m | 119 +++++++++++++++++++++++++++++++++ +rsa/+fsl/noiseNormalizeBeta.m | 109 ++++++++++++++++++++++++++++++ 2 files changed, 228 insertions(+) create mode 100644 +rsa/+fsl/getDataFromFSL.m create mode 100644 +rsa/+fsl/noiseNormalizeBeta.m diff --git a/+rsa/+fsl/getDataFromFSL.m b/+rsa/+fsl/getDataFromFSL.m new file mode 100644 index 0000000..ce0fde8 --- /dev/null +++ b/+rsa/+fsl/getDataFromFSL.m @@ -0,0 +1,119 @@ +function [betas,residuals] = getDataFromFSL(userOptions) +% +% getDataFromSPM is a function which will extract from the SPM metadata the +% correspondence between the beta image filenames and the condition and session +% number. +% +% function betas = getDataFromSPM(userOptions) +% +% betas --- The array of info. +% betas(condition, session).identifier is a string which referrs +% to the filename (not including path) of the SPM beta image. +% +% userOptions --- The options struct. +% userOptions.subjectNames +% A cell array containing strings identifying the subject +% names. +% userOptions.betaPath +% A string which contains the absolute path to the +% location of the beta images. It can contain the +% following wildcards which would be replaced as +% indicated: +% [[subjectName]] +% To be replaced with the contents of +% subject. +% [[betaIdentifier]] +% To be replaced by filenames returned in +% betas. +% userOptions.conditionLabels +% A cell array containing the names of the conditions in +% this experiment. Here, these will be used to find +% condition response beta predictor images (so make sure +% they're the same as the ones used for SPM!). +% +% Cai Wingfield 12-2009, 6-2010, 8-2010 + + import rsa.* +import rsa.fig.* +import rsa.fmri.* +import rsa.rdm.* +import rsa.sim.* +import rsa.spm.* +import rsa.stat.* +import rsa.util.* + + %% Set defaults and check for problems. + if ~isfield(userOptions, 'betaPath'), error('getDataFromSPM:NoBetaPath', 'userOptions.betaPath is not set. See help.'); end%if + if ~isfield(userOptions, 'subjectNames'), error('getDataFromSPM:NoSubjectNames', 'userOptions.subjectNames is not set. See help.'); end%if + if ~isfield(userOptions, 'conditionLabels'), error('getDataFromSPM:NoConditionLabels', 'userOptions.conditionLabels is not set. See help.'); end%if + + firstSubject = userOptions.subjectNames{1}; + + readFile = replaceWildcards(fullfile(userOptions.betaPath, 'SPM.mat'), '[[subjectName]]', firstSubject, '[[betaIdentifier]]', ''); + load(readFile); + nBetas = max(size(SPM.Vbeta)); + + nConditions = numel(userOptions.conditionLabels); + + % Extract all info + + highestSessionNumber = 0; + + for betaNumber = 1:nBetas % For each beta file... + + % Get the description of the beta file + thisBetaDescrip = SPM.Vbeta(betaNumber).descrip; + + % Extract from it the file name, the session number and the condition name + [thisBetaName, thisSessionNumber, thisConditionName] = extractSingleBetaInfo(thisBetaDescrip); + + % Check if it's one of the conditions + thisBetaIsACondition = false; % (What a delicious variable name!) + for condition = 1:nConditions + + thisBetaIsACondition = strcmpi(thisConditionName, userOptions.conditionLabels{condition}); % Check if it's a condition specified by the user + conditionNumber = condition; % Record the condition number that it is + if thisBetaIsACondition, break; end%if + + end%for + + % Now only proceed if we're looking at a condition + if thisBetaIsACondition + + % Keep a tally of the higest run yet + highestSessionNumber = max(highestSessionNumber, thisSessionNumber); + + % Store the file name in the betas struct + betas(thisSessionNumber, conditionNumber).identifier = [thisBetaName '.img']; + + end%if + end%for + +end%function + +%% Subfunctions: %% + +% spm_spm:beta (0001) - Sn(1) all_events*bf(1) +function [betaName, sessionNumber, conditionName] = extractSingleBetaInfo(strIn) + + import rsa.* +import rsa.fig.* +import rsa.fmri.* +import rsa.rdm.* +import rsa.sim.* +import rsa.spm.* +import rsa.stat.* +import rsa.util.* + + openBrackets = findstr(strIn, '('); + closedBrackets = findstr(strIn, ')'); + star = findstr(strIn, '*'); + + inFirstBrackets = strIn(openBrackets(1)+1:closedBrackets(1)-1); + betaName = ['beta_' inFirstBrackets]; + + inSecondBrackets = strIn(openBrackets(2)+1:closedBrackets(2)-1); + sessionNumber = str2num(inSecondBrackets); + + conditionName = strIn(closedBrackets(2)+2:star-1); +end%function diff --git a/+rsa/+fsl/noiseNormalizeBeta.m b/+rsa/+fsl/noiseNormalizeBeta.m new file mode 100644 index 0000000..b420fae --- /dev/null +++ b/+rsa/+fsl/noiseNormalizeBeta.m @@ -0,0 +1,109 @@ +function [u_hat,resMS,Sw_hat,beta_hat,shrinkage,trRR]=noiseNormalizeBeta(Y,SPM,varargin) +% function [u_hat,Sw_hat,resMS,beta_hat]=rsa_noiseNormalizeBeta(Y,SPM,varargin) +% Estimates beta coefficiencts beta_hat and residuals from raw time series Y +% Estimates the true activity patterns u_hat by applying noise normalization to beta_hat +% INPUT: +% Y raw timeseries, T by P +% SPM: SPM structure +% OPTIONS: +% 'normmode': 'overall': Does the multivariate noise normalisation overall (default) +% 'runwise': Does the multivariate noise normalisation by run +% 'shrinkage': Shrinkage coefficient. +% 0: No regularisation +% 1: Using only the diagonal - i.e. univariate noise normalisation +% By default the shrinkage coeffcient is determined using +% the Ledoit-Wolf method. +% OUTPUT: +% u_hat: estimated true activity patterns (beta_hat after multivariate noise normalization), +% resMS: residual mean-square - diagonal of the Var-cov matrix, 1 by P +% Sw_hat: overall voxel error variance-covariance matrix (PxP), before regularisation +% beta_hat: estimated raw regression coefficients, K*R by P +% shrinkage: applied shrinkage factor +% trRR: Trace of the squared residual spatial correlation matrix +% This number provides an estimate for the residual spatial +% correlation. For variance-estimation, the effective +% number of voxels are effVox = numVox^2/trRR +% Alexander Walther, Joern Diedrichsen +% joern.diedrichsen@googlemail.com +% 2/2015 +Opt.normmode = 'overall'; % Either runwise or overall +Opt.shrinkage = []; +Opt = rsa.getUserOptions(varargin,Opt); +[T,P]=size(Y); %%% number of time points and voxels + +%%% Discard NaN voxels +test=isnan(sum(Y)); +if (any(test)) + warning(sprintf('%d of %d voxels contained NaNs -discarding',sum(test),length(test))); + Y=Y(:,test==0); +end; + +xX = SPM.xX; %%% take the design +[T,Q] = size(xX.X); + +%%% Get partions: For each run (1:K), find the time points (T) and regressors (K+Q) that belong to the run +partT = nan(T,1); +partQ = nan(Q,1); +Nrun=length(SPM.Sess); %%% number of runs +for i=1:Nrun + partT(SPM.Sess(i).row,1)=i; + partQ(SPM.Sess(i).col,1)=i; + partQ(SPM.xX.iB(i),1)=i; %%% Add intercepts +end; + +%%% redo the first-level GLM using matlab functions +KWY=spm_filter(xX.K,xX.W*Y); %%% filter out low-frequence trends in Y +beta_hat=xX.pKX*KWY; %%% ordinary least squares estimate of beta_hat = inv(X'*X)*X'*Y +res=spm_sp('r',xX.xKXs,KWY); %%% residuals: res = Y - X*beta +clear KWY XZ %%% clear to save memory + +switch (Opt.normmode) + case 'runwise' % do run-wise noise normalization + u_hat = zeros(size(beta_hat)); + shrink=zeros(Nrun,1); + for i=1:Nrun + idxT = partT==i; % Time points for this partition + idxQ = partQ==i; % Regressors for this partition + numFilt = size(xX.K(i).X0,2); % Number of filter variables for this run + + [Sw_reg(:,:,i),shrinkage(i),Sw_hat(:,:,i)]=rsa.stat.covdiag(res(idxT,:),sum(idxT)-sum(idxQ)-numFilt-1,'shrinkage',Opt.shrinkage); %%% regularize Sw_hat through optimal shrinkage + % Postmultiply by the inverse square root of the estimated matrix + [V,L]=eig(Sw_reg(:,:,i)); % This is overall faster and numerical more stable than Sw_hat.^-1/2 + l=diag(L); + sq(:,:,i) = V*bsxfun(@rdivide,V',sqrt(l)); % Slightly faster than sq = V*diag(1./sqrt(l))*V'; + u_hat(idxQ,:)=beta_hat(idxQ,:)*sq(:,:,i); + end; + shrinkage=mean(shrinkage); + Sw_hat = mean(Sw_hat,3); + Sw_reg = mean(Sw_reg,3); + case 'overall' %%% do overall noise normalization + [Sw_reg,shrinkage,Sw_hat]=rsa.stat.covdiag(res,SPM.xX.erdf,'shrinkage',Opt.shrinkage); %%% regularize Sw_hat through optimal shrinkage + + % Postmultiply by the inverse square root of the estimated matrix + [V,L]=eig(Sw_reg); + l=diag(L); + sq = V*bsxfun(@rdivide,V',sqrt(l)); % Slightly faster than sq = V*diag(1./sqrt(l))*V'; + u_hat=beta_hat*sq; +end; + +% Return the diagonal of Sw_hat +if (nargout>1) + resMS=sum(res.^2)./SPM.xX.erdf; +end; + +% return the trace of the residual covariance matrix +% If prewhiting was perfect, then the trace(sig*sig) = P and the effective number of voxels +% would be P. However, because we need to regularise our estimate, the residual +% spatial covariance matrix is somwhat correlated. For variance estimation on the distances, +% the effective voxel number would be numVox^2/trRR +if (nargout>5) + % Caluclate the predicted redidual covariance + sq = mean(sq,3); + V_hat = sq'*Sw_hat*sq; + % Calculate the correlation matrix from this + Vsigma = sqrt(diag(V_hat)); + R_hat = bsxfun(@rdivide,V_hat,Vsigma); + R_hat = bsxfun(@rdivide,R_hat,Vsigma'); % R = C ./ sigma*sigma'; + % get the trace + trRR = sum(sum(R_hat.*R_hat)); % traceRR +end; From a2b6e56c2e5e1b9c6e6c8e1db11fa7705decf4be Mon Sep 17 00:00:00 2001 From: Roni Maimon Date: Wed, 6 Sep 2017 15:35:47 +0100 Subject: [PATCH 3/7] Added the config file for the fsl demo --- Demos/projectOptions_fsl_demo.m | 74 +++++++++++++++++++++++++++++++++ 1 file changed, 74 insertions(+) create mode 100644 Demos/projectOptions_fsl_demo.m diff --git a/Demos/projectOptions_fsl_demo.m b/Demos/projectOptions_fsl_demo.m new file mode 100644 index 0000000..ae801c2 --- /dev/null +++ b/Demos/projectOptions_fsl_demo.m @@ -0,0 +1,74 @@ +function userOptions = projectOptions_fsl_demo() +% projectOptions_fsl_demo is an options file for the FSL demo tutorial +% +% Roni Maimon 9-2017 +%__________________________________________________________________________ +% Copyright (C) + +%% Project details +% This name identifies a collection of files which all belong to the same run of a project. +userOptions.analysisName = 'DEMO_FSL'; % this is renamed in the code for demos 3-4. + +% This is the root directory of the project. +userOptions.rootPath = [pwd,filesep,'DEMO_FSL']; + + +%%%%%%%%%%%%%%%%%%%%%%%% +%% EXPERIMENTAL SETUP %% +%%%%%%%%%%%%%%%%%%%%%%%% + +% The list of subjects to be included in the study. +userOptions.subjectNames = { ... + 'T05', ... + 'T06', ... + }; + +%% FSL Specific parameters +userOptions.run_names = { ... + 'A1b','A2b','A3b','A4b' ... + }; +userOptions.featsPrefix = '5_'; +userOptions.featsSuffix = ''; +userOptions.featsPath = [pwd,filesep,'FSLData',filesep,'[[subjectName]]',... + filesep,'models',filesep,'glm',filesep,'fingers', ... + filesep, '[[featPrefix]]','[[runName]]','[[featSufix]]','.feat']; + +userOptions.copes = { ... + 1,2,3,4,5 ... + }; + +%% End of FSL specific parameters + +%%%%%%%%%%%%%%%%%%%%%%%%%% +%% ANALYSIS PREFERENCES %% +%%%%%%%%%%%%%%%%%%%%%%%%%% + +%% First-order analysis +% Text lables which may be attached to the conditions for MDS plots. +userOptions.conditionLabels = { ... + 'F1', ... + 'F2', ... + 'F3', ... + 'F4', ... + 'F5' + }; + +userOptions.useAlternativeConditionLabels = false; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% FEATUERS OF INTEREST SELECTION OPTIONS %% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + %% %% %% %% %% + %% fMRI %% Use these next three options if you're working in fMRI native space: + %% %% %% %% %% + + % The path to a stereotypical mask data file is stored (not including subject-specific identifiers). + % "[[subjectName]]" should be used as a placeholder to denote an entry in userOptions.subjectNames + % "[[maskName]]" should be used as a placeholder to denote an entry in userOptions.maskNames + userOptions.maskPath = [pwd,filesep,'FSLData',filesep,'[[subjectName]]',filesep,'models',filesep,'rsa',filesep,'masks',filesep,'[[maskName]].nii.gz']; + + % The list of mask filenames (minus .hdr extension) to be used. + userOptions.maskNames = { ... + 'Left_S1_handknob_alldigits'... + }; +end%function From 7aaa46b693cd610623115f70065b4974cf5e677b Mon Sep 17 00:00:00 2001 From: Roni Maimon Date: Thu, 7 Sep 2017 18:51:46 +0100 Subject: [PATCH 4/7] Added FSL demo data link --- Demos/FSLData/Readme.md | 3 +++ Demos/projectOptions_fsl_demo.m | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) create mode 100644 Demos/FSLData/Readme.md diff --git a/Demos/FSLData/Readme.md b/Demos/FSLData/Readme.md new file mode 100644 index 0000000..c82d3b4 --- /dev/null +++ b/Demos/FSLData/Readme.md @@ -0,0 +1,3 @@ +# Download instructions +Downlad FSLData From:https://www.dropbox.com/s/usjjz31tf3q5kwx/FSLData.zip?dl=0 +and extract here \ No newline at end of file diff --git a/Demos/projectOptions_fsl_demo.m b/Demos/projectOptions_fsl_demo.m index ae801c2..fb5e472 100644 --- a/Demos/projectOptions_fsl_demo.m +++ b/Demos/projectOptions_fsl_demo.m @@ -69,6 +69,6 @@ % The list of mask filenames (minus .hdr extension) to be used. userOptions.maskNames = { ... - 'Left_S1_handknob_alldigits'... + 'S1_handknob_alldigits'... }; end%function From 87226324a801d23e44ef272e0990cdbda86014af Mon Sep 17 00:00:00 2001 From: Daan Wesselink Date: Fri, 17 Nov 2017 14:35:29 +0000 Subject: [PATCH 5/7] loading fsl data and creating distance matrix for single subject --- +rsa/+fsl/getDataFromFSL.m | 151 +++++++----------------------- +rsa/+fsl/noiseNormalizeBetaFSL.m | 17 ++++ Demos/defineUserOptions.m | 2 +- Demos/fsl_demo.m | 17 ++++ Demos/projectOptions_fsl_demo.m | 12 ++- 5 files changed, 75 insertions(+), 124 deletions(-) create mode 100644 +rsa/+fsl/noiseNormalizeBetaFSL.m create mode 100644 Demos/fsl_demo.m diff --git a/+rsa/+fsl/getDataFromFSL.m b/+rsa/+fsl/getDataFromFSL.m index ce0fde8..62bae32 100644 --- a/+rsa/+fsl/getDataFromFSL.m +++ b/+rsa/+fsl/getDataFromFSL.m @@ -1,119 +1,32 @@ -function [betas,residuals] = getDataFromFSL(userOptions) -% -% getDataFromSPM is a function which will extract from the SPM metadata the -% correspondence between the beta image filenames and the condition and session -% number. -% -% function betas = getDataFromSPM(userOptions) -% -% betas --- The array of info. -% betas(condition, session).identifier is a string which referrs -% to the filename (not including path) of the SPM beta image. -% -% userOptions --- The options struct. -% userOptions.subjectNames -% A cell array containing strings identifying the subject -% names. -% userOptions.betaPath -% A string which contains the absolute path to the -% location of the beta images. It can contain the -% following wildcards which would be replaced as -% indicated: -% [[subjectName]] -% To be replaced with the contents of -% subject. -% [[betaIdentifier]] -% To be replaced by filenames returned in -% betas. -% userOptions.conditionLabels -% A cell array containing the names of the conditions in -% this experiment. Here, these will be used to find -% condition response beta predictor images (so make sure -% they're the same as the ones used for SPM!). -% -% Cai Wingfield 12-2009, 6-2010, 8-2010 - - import rsa.* -import rsa.fig.* -import rsa.fmri.* -import rsa.rdm.* -import rsa.sim.* -import rsa.spm.* -import rsa.stat.* -import rsa.util.* - - %% Set defaults and check for problems. - if ~isfield(userOptions, 'betaPath'), error('getDataFromSPM:NoBetaPath', 'userOptions.betaPath is not set. See help.'); end%if - if ~isfield(userOptions, 'subjectNames'), error('getDataFromSPM:NoSubjectNames', 'userOptions.subjectNames is not set. See help.'); end%if - if ~isfield(userOptions, 'conditionLabels'), error('getDataFromSPM:NoConditionLabels', 'userOptions.conditionLabels is not set. See help.'); end%if - - firstSubject = userOptions.subjectNames{1}; - - readFile = replaceWildcards(fullfile(userOptions.betaPath, 'SPM.mat'), '[[subjectName]]', firstSubject, '[[betaIdentifier]]', ''); - load(readFile); - nBetas = max(size(SPM.Vbeta)); - - nConditions = numel(userOptions.conditionLabels); - - % Extract all info - - highestSessionNumber = 0; - - for betaNumber = 1:nBetas % For each beta file... - - % Get the description of the beta file - thisBetaDescrip = SPM.Vbeta(betaNumber).descrip; - - % Extract from it the file name, the session number and the condition name - [thisBetaName, thisSessionNumber, thisConditionName] = extractSingleBetaInfo(thisBetaDescrip); - - % Check if it's one of the conditions - thisBetaIsACondition = false; % (What a delicious variable name!) - for condition = 1:nConditions - - thisBetaIsACondition = strcmpi(thisConditionName, userOptions.conditionLabels{condition}); % Check if it's a condition specified by the user - conditionNumber = condition; % Record the condition number that it is - if thisBetaIsACondition, break; end%if - - end%for - - % Now only proceed if we're looking at a condition - if thisBetaIsACondition - - % Keep a tally of the higest run yet - highestSessionNumber = max(highestSessionNumber, thisSessionNumber); - - % Store the file name in the betas struct - betas(thisSessionNumber, conditionNumber).identifier = [thisBetaName '.img']; - - end%if - end%for - -end%function - -%% Subfunctions: %% - -% spm_spm:beta (0001) - Sn(1) all_events*bf(1) -function [betaName, sessionNumber, conditionName] = extractSingleBetaInfo(strIn) - - import rsa.* -import rsa.fig.* -import rsa.fmri.* -import rsa.rdm.* -import rsa.sim.* -import rsa.spm.* -import rsa.stat.* -import rsa.util.* - - openBrackets = findstr(strIn, '('); - closedBrackets = findstr(strIn, ')'); - star = findstr(strIn, '*'); - - inFirstBrackets = strIn(openBrackets(1)+1:closedBrackets(1)-1); - betaName = ['beta_' inFirstBrackets]; - - inSecondBrackets = strIn(openBrackets(2)+1:closedBrackets(2)-1); - sessionNumber = str2num(inSecondBrackets); - - conditionName = strIn(closedBrackets(2)+2:star-1); -end%function +function [betas,residuals] = getDataFromFSL(userOptions,subjectName) + +%load mask +curmaskPath=replaceWildcards(userOptions.maskPath, ... + '[[subjectName]]',subjectName, ... + '[[maskName]]',userOptions.maskNames{1}); +mask=load_nii(curmaskPath); + +%load copes +betas=[]; +residuals=[]; +for r = userOptions.run_names + curfeatsPath=replaceWildcards(userOptions.featsPath, ... + '[[subjectName]]',subjectName, ... + '[[featPrefix]]',userOptions.featsPrefix, ... + '[[runName]]',r{1}, ... + '[[featSuffix]]',userOptions.featsSuffix); + for c = userOptions.copes + beta=load_nii(fullfile(curfeatsPath,'stats',['cope' int2str(c{1})' '.nii.gz'])); + betas=cat(2,betas,beta.img(boolean(mask.img))); + end + + res=load_nii(fullfile(curfeatsPath,'stats','res4d.nii.gz')); + res1=res.img(repmat(boolean(mask.img),[1 1 1 size(res.img,4)])); + residuals=cat(2,residuals,reshape(res1,[length(res1)/size(res.img,4) size(res.img,4)])); + %ress=cat(2,ress,res.img(boolean(mask.img))); +end + +residuals=residuals'; +betas=betas'; + +end diff --git a/+rsa/+fsl/noiseNormalizeBetaFSL.m b/+rsa/+fsl/noiseNormalizeBetaFSL.m new file mode 100644 index 0000000..d790bf0 --- /dev/null +++ b/+rsa/+fsl/noiseNormalizeBetaFSL.m @@ -0,0 +1,17 @@ +function [u_hat] = noiseNormalizeBetaFSL(Beta,Res,Partition) + +[T,P] = size(Res); +[Q,P] = size(Beta); + +Nrun= max(Partition); +partT = kron(1:Nrun,ones(1,T/Nrun))'; +partQ = Partition; + +for i=1:Nrun + idxT=partT==i; + idxQ=partQ==i; + Sw_hat(:,:,i)=rsa.stat.covdiag(Res(idxT,:)); %%% regularize Sw_hat through optimal shrinkage + u_hat(idxQ,:)=Beta(idxQ,:)*Sw_hat(:,:,i)^(-1/2); %%% multivariate noise normalization +end; + +end \ No newline at end of file diff --git a/Demos/defineUserOptions.m b/Demos/defineUserOptions.m index 1882876..5f66448 100644 --- a/Demos/defineUserOptions.m +++ b/Demos/defineUserOptions.m @@ -8,7 +8,7 @@ % a new name, so all will not be lost if you don't do this). % % For a guide to how to fill out the fields in this file, consult -% the documentation folder (particularly the userOptions_guide.m) +% the docu mentation folder (particularly the userOptions_guide.m) % % Cai Wingfield 11-2009 %__________________________________________________________________________ diff --git a/Demos/fsl_demo.m b/Demos/fsl_demo.m new file mode 100644 index 0000000..89449e1 --- /dev/null +++ b/Demos/fsl_demo.m @@ -0,0 +1,17 @@ +function distance=fsl_demo() + +userOptions=projectOptions_fsl_demo; + +[betas,residuals]=rsa.fsl.getDataFromFSL(userOptions,'T05'); + +partQ = kron(1:length(userOptions.run_names),ones(1,length(userOptions.copes)))'; %runNumber (partition) +partK = repmat(1:length(userOptions.copes),[1 length(userOptions.run_names)])'; %condNumber (CondVec) + +u_hat = rsa.fsl.noiseNormalizeBetaFSL(betas,residuals,partQ); % Get noise normalised betas +distance = rsa.distanceLDC(u_hat,partQ,partK); + +figure; imagesc(squareform(distance)); + + + +end \ No newline at end of file diff --git a/Demos/projectOptions_fsl_demo.m b/Demos/projectOptions_fsl_demo.m index fb5e472..c34a3b7 100644 --- a/Demos/projectOptions_fsl_demo.m +++ b/Demos/projectOptions_fsl_demo.m @@ -6,17 +6,21 @@ % Copyright (C) %% Project details +%Current folder +curFolder=fileparts(which('projectOptions_fsl_demo.m')); + % This name identifies a collection of files which all belong to the same run of a project. userOptions.analysisName = 'DEMO_FSL'; % this is renamed in the code for demos 3-4. % This is the root directory of the project. -userOptions.rootPath = [pwd,filesep,'DEMO_FSL']; +userOptions.rootPath = [curFolder,filesep,'DEMO_FSL']; %%%%%%%%%%%%%%%%%%%%%%%% %% EXPERIMENTAL SETUP %% %%%%%%%%%%%%%%%%%%%%%%%% + % The list of subjects to be included in the study. userOptions.subjectNames = { ... 'T05', ... @@ -29,9 +33,9 @@ }; userOptions.featsPrefix = '5_'; userOptions.featsSuffix = ''; -userOptions.featsPath = [pwd,filesep,'FSLData',filesep,'[[subjectName]]',... +userOptions.featsPath = [curFolder,filesep,'FSLData',filesep,'[[subjectName]]',... filesep,'models',filesep,'glm',filesep,'fingers', ... - filesep, '[[featPrefix]]','[[runName]]','[[featSufix]]','.feat']; + filesep, '[[featPrefix]]','[[runName]]','[[featSuffix]]','.feat']; userOptions.copes = { ... 1,2,3,4,5 ... @@ -65,7 +69,7 @@ % The path to a stereotypical mask data file is stored (not including subject-specific identifiers). % "[[subjectName]]" should be used as a placeholder to denote an entry in userOptions.subjectNames % "[[maskName]]" should be used as a placeholder to denote an entry in userOptions.maskNames - userOptions.maskPath = [pwd,filesep,'FSLData',filesep,'[[subjectName]]',filesep,'models',filesep,'rsa',filesep,'masks',filesep,'[[maskName]].nii.gz']; + userOptions.maskPath = [curFolder,filesep,'FSLData',filesep,'[[subjectName]]',filesep,'models',filesep,'masks',filesep,'[[maskName]].nii.gz']; % The list of mask filenames (minus .hdr extension) to be used. userOptions.maskNames = { ... From 729a351622029d7d42dd78b3e85c6ce647121a76 Mon Sep 17 00:00:00 2001 From: Daan Wesselink Date: Fri, 17 Nov 2017 14:38:45 +0000 Subject: [PATCH 6/7] using right replaceWildcards --- +rsa/+fsl/getDataFromFSL.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/+rsa/+fsl/getDataFromFSL.m b/+rsa/+fsl/getDataFromFSL.m index 62bae32..b01a4e5 100644 --- a/+rsa/+fsl/getDataFromFSL.m +++ b/+rsa/+fsl/getDataFromFSL.m @@ -1,7 +1,7 @@ function [betas,residuals] = getDataFromFSL(userOptions,subjectName) %load mask -curmaskPath=replaceWildcards(userOptions.maskPath, ... +curmaskPath=rsa.util.replaceWildcards(userOptions.maskPath, ... '[[subjectName]]',subjectName, ... '[[maskName]]',userOptions.maskNames{1}); mask=load_nii(curmaskPath); @@ -10,7 +10,7 @@ betas=[]; residuals=[]; for r = userOptions.run_names - curfeatsPath=replaceWildcards(userOptions.featsPath, ... + curfeatsPath=rsa.util.replaceWildcards(userOptions.featsPath, ... '[[subjectName]]',subjectName, ... '[[featPrefix]]',userOptions.featsPrefix, ... '[[runName]]',r{1}, ... From 20fbe053501a1fa891cd16f60a8c07b6ec12d168 Mon Sep 17 00:00:00 2001 From: Roni Maimon Date: Sun, 26 Nov 2017 15:25:54 +0000 Subject: [PATCH 7/7] Removing unnecessary files --- +rsa/+fsl/noiseNormalizeBeta.m | 109 --------------------------------- Hello/test.m | 1 - 2 files changed, 110 deletions(-) delete mode 100644 +rsa/+fsl/noiseNormalizeBeta.m delete mode 100644 Hello/test.m diff --git a/+rsa/+fsl/noiseNormalizeBeta.m b/+rsa/+fsl/noiseNormalizeBeta.m deleted file mode 100644 index b420fae..0000000 --- a/+rsa/+fsl/noiseNormalizeBeta.m +++ /dev/null @@ -1,109 +0,0 @@ -function [u_hat,resMS,Sw_hat,beta_hat,shrinkage,trRR]=noiseNormalizeBeta(Y,SPM,varargin) -% function [u_hat,Sw_hat,resMS,beta_hat]=rsa_noiseNormalizeBeta(Y,SPM,varargin) -% Estimates beta coefficiencts beta_hat and residuals from raw time series Y -% Estimates the true activity patterns u_hat by applying noise normalization to beta_hat -% INPUT: -% Y raw timeseries, T by P -% SPM: SPM structure -% OPTIONS: -% 'normmode': 'overall': Does the multivariate noise normalisation overall (default) -% 'runwise': Does the multivariate noise normalisation by run -% 'shrinkage': Shrinkage coefficient. -% 0: No regularisation -% 1: Using only the diagonal - i.e. univariate noise normalisation -% By default the shrinkage coeffcient is determined using -% the Ledoit-Wolf method. -% OUTPUT: -% u_hat: estimated true activity patterns (beta_hat after multivariate noise normalization), -% resMS: residual mean-square - diagonal of the Var-cov matrix, 1 by P -% Sw_hat: overall voxel error variance-covariance matrix (PxP), before regularisation -% beta_hat: estimated raw regression coefficients, K*R by P -% shrinkage: applied shrinkage factor -% trRR: Trace of the squared residual spatial correlation matrix -% This number provides an estimate for the residual spatial -% correlation. For variance-estimation, the effective -% number of voxels are effVox = numVox^2/trRR -% Alexander Walther, Joern Diedrichsen -% joern.diedrichsen@googlemail.com -% 2/2015 -Opt.normmode = 'overall'; % Either runwise or overall -Opt.shrinkage = []; -Opt = rsa.getUserOptions(varargin,Opt); -[T,P]=size(Y); %%% number of time points and voxels - -%%% Discard NaN voxels -test=isnan(sum(Y)); -if (any(test)) - warning(sprintf('%d of %d voxels contained NaNs -discarding',sum(test),length(test))); - Y=Y(:,test==0); -end; - -xX = SPM.xX; %%% take the design -[T,Q] = size(xX.X); - -%%% Get partions: For each run (1:K), find the time points (T) and regressors (K+Q) that belong to the run -partT = nan(T,1); -partQ = nan(Q,1); -Nrun=length(SPM.Sess); %%% number of runs -for i=1:Nrun - partT(SPM.Sess(i).row,1)=i; - partQ(SPM.Sess(i).col,1)=i; - partQ(SPM.xX.iB(i),1)=i; %%% Add intercepts -end; - -%%% redo the first-level GLM using matlab functions -KWY=spm_filter(xX.K,xX.W*Y); %%% filter out low-frequence trends in Y -beta_hat=xX.pKX*KWY; %%% ordinary least squares estimate of beta_hat = inv(X'*X)*X'*Y -res=spm_sp('r',xX.xKXs,KWY); %%% residuals: res = Y - X*beta -clear KWY XZ %%% clear to save memory - -switch (Opt.normmode) - case 'runwise' % do run-wise noise normalization - u_hat = zeros(size(beta_hat)); - shrink=zeros(Nrun,1); - for i=1:Nrun - idxT = partT==i; % Time points for this partition - idxQ = partQ==i; % Regressors for this partition - numFilt = size(xX.K(i).X0,2); % Number of filter variables for this run - - [Sw_reg(:,:,i),shrinkage(i),Sw_hat(:,:,i)]=rsa.stat.covdiag(res(idxT,:),sum(idxT)-sum(idxQ)-numFilt-1,'shrinkage',Opt.shrinkage); %%% regularize Sw_hat through optimal shrinkage - % Postmultiply by the inverse square root of the estimated matrix - [V,L]=eig(Sw_reg(:,:,i)); % This is overall faster and numerical more stable than Sw_hat.^-1/2 - l=diag(L); - sq(:,:,i) = V*bsxfun(@rdivide,V',sqrt(l)); % Slightly faster than sq = V*diag(1./sqrt(l))*V'; - u_hat(idxQ,:)=beta_hat(idxQ,:)*sq(:,:,i); - end; - shrinkage=mean(shrinkage); - Sw_hat = mean(Sw_hat,3); - Sw_reg = mean(Sw_reg,3); - case 'overall' %%% do overall noise normalization - [Sw_reg,shrinkage,Sw_hat]=rsa.stat.covdiag(res,SPM.xX.erdf,'shrinkage',Opt.shrinkage); %%% regularize Sw_hat through optimal shrinkage - - % Postmultiply by the inverse square root of the estimated matrix - [V,L]=eig(Sw_reg); - l=diag(L); - sq = V*bsxfun(@rdivide,V',sqrt(l)); % Slightly faster than sq = V*diag(1./sqrt(l))*V'; - u_hat=beta_hat*sq; -end; - -% Return the diagonal of Sw_hat -if (nargout>1) - resMS=sum(res.^2)./SPM.xX.erdf; -end; - -% return the trace of the residual covariance matrix -% If prewhiting was perfect, then the trace(sig*sig) = P and the effective number of voxels -% would be P. However, because we need to regularise our estimate, the residual -% spatial covariance matrix is somwhat correlated. For variance estimation on the distances, -% the effective voxel number would be numVox^2/trRR -if (nargout>5) - % Caluclate the predicted redidual covariance - sq = mean(sq,3); - V_hat = sq'*Sw_hat*sq; - % Calculate the correlation matrix from this - Vsigma = sqrt(diag(V_hat)); - R_hat = bsxfun(@rdivide,V_hat,Vsigma); - R_hat = bsxfun(@rdivide,R_hat,Vsigma'); % R = C ./ sigma*sigma'; - % get the trace - trRR = sum(sum(R_hat.*R_hat)); % traceRR -end; diff --git a/Hello/test.m b/Hello/test.m deleted file mode 100644 index 5ab2f8a..0000000 --- a/Hello/test.m +++ /dev/null @@ -1 +0,0 @@ -Hello \ No newline at end of file