This should be the alpha version.

resources/BSM_Class/BSM.m

@@ -0,0 +1,451 @@
+classdef BSM<handle
+    %BSM Summary of this class goes here
+    %   Detailed explanation goes here
+    
+    properties %public properties
+        SPMs={}; %field to contain the full path of all SPM.mat files to be analysed
+        Acfl %columns in SPM.xX.X to ignore
+        Type=4; %BSM type 1..4 5: none %default is option 4
+        KeepMP = false; %retain Motion paramters
+        Prefix = 'bs'; %beta series prefix
+        Clean = true;
+    end
+    properties (Access = private)%private properties
+        SPM_orig %temporary storage for SPM.mat file
+        ActiveColumns %store the columns to be analized
+        U_orig %store the U field of SPM
+        TargetFolder %target folder
+        % option below obsolete?
+        BSMToolboxPath%store path to BSMToolbox
+        SPMPath;%store path to SPM folder
+        SPMIndexList=[];
+        MonitorState=false;
+    end
+    properties (Dependent) %dependent properties
+        NumSPMs% number of SPM.mat files to process
+        
+    end
+    properties (Constant = true)%constants
+        BSMMethod={'Traditional BSM Rissman et al', 'Adapted BSM Rissman et al', 'Fast approach (Mumford et al)', 'Adapted fast approach','none'};
+    end
+    
+    methods %user callable methods
+        function Monitor(obj,varargin)
+            if numel(varargin)>1
+                id = 'BSM:Monitor:TooManyInputs';
+                msg= sprintf('Too many inputs given.\n Use: \t obj.Monitor\n\tobj.Monitor(true)\n\tobj.monitor(false)\n');
+                throw(MException(id,msg))
+            end
+            if numel(varargin)==0
+                obj.MonitorState=xor(obj.MonitorState,true);
+            else
+                if ~islogical(varargin{1})
+                    id = 'BSM:Monitor:InputIsNotLogical';
+                    msg= sprintf('Input is not of type logical\n');
+                    throw(MException(id,msg))
+                end
+                obj.MonitorState=varargin{1};
+            end
+        end
+        function Calc(obj,varargin)
+            % obj.calc to apply calculate to all SPM.mat files
+            % obj.calc(n) to calculated the nth SPM.mat file
+            % obj.calc([1 3 6]) to indicate which SPM.mat file to calculate
+            
+            %only initialize once. Every subsequent call to this function,
+            %do not reinitialize. 
+            persistent InitDone
+            if isempty(InitDone) || not(InitDone)
+                %   initialize spm job enviourment
+                spm_jobman('initcfg');
+                InitDone = true;
+            end
+            SPMidx=[];% will contain the index to the SPMs to work on
+            if logical(numel(varargin))
+                SPMidx=varargin{1};
+            end
+            %check number of input arguments
+            if ~logical(obj.NumSPMs)%nothing to do
+                fprintf('== == == Nothing to do: done == == ==\n')
+                return;
+            end
+            if isempty(SPMidx)
+                SPMidx=1:obj.NumSPMs;
+            end
+            if ~isnumeric(SPMidx)||~all(isfinite(SPMidx))||min(SPMidx)<1||max(SPMidx)>obj.NumSPMs
+                id='BSM:Calc:InvalidIndex';
+                msg='Invalid Index';
+                throw(MException(id,msg));
+            end
+            if numel(SPMidx)>1
+                obj.SPMIndexList=SPMidx;
+                SPMidx=obj.SPMIndexList(end);
+                try
+                    obj.Calc(SPMidx); %handle the single case
+                catch ME
+                    %if the identifier matches, just continue with next SPM
+                    if ~any(strcmpi(ME.identifier,...
+                            {'BSM:CoreCalc:SPMxBFOrder:invalid'...
+                            ,'BSM:CoreCalc:SPMSess:invalid'...
+                            }))
+                        rethrow(ME)
+                    else
+                        disp(ME.message);
+                    end
+                end
+                obj.SPMIndexList=obj.SPMIndexList(1:end-1);
+                if isempty(obj.SPMIndexList)
+                    return;
+                else
+                    obj.Calc(obj.SPMIndexList);% recursively run the rest of the list
+                end
+            else %  deal with single entry
+                obj.CoreCalc(SPMidx)
+            end
+            
+            
+        end
+    end
+    methods %constructor and set/get methods
+        function obj = BSM(varargin)
+            %BSM Construct an instance of this class
+            %option: parse the info of a matlab batch file?  see
+            %tbx_cf_bsm?
+            %TODO:
+            % set default values
+            obj.SPMPath=[]; %get the SPM path
+            if nargin>0
+                [val,TF]=obj.MyGetVal('SPMs',varargin);if TF;obj.SPMs=val;end
+                [val,TF]=obj.MyGetVal('Type',varargin);if TF;obj.Type=val;end
+                [val,TF]=obj.MyGetVal('Acfl',varargin);if TF;obj.Acfl=val;end
+            end
+           
+        end
+        function val=get.NumSPMs(obj)
+            val=numel(obj.SPMs);
+        end
+        function set.SPMPath(obj,~)
+            SPMFound=which('spm'); %check if SPM in searchpath
+            if isempty(SPMFound) %SPM not found
+                error('please add SPM to your path before trying to run BSM');
+                %to do add a check on spm version
+            end
+            obj.SPMPath=fileparts(SPMFound);
+        end
+        function set.Clean(obj,varargin)
+            if isempty(varargin)
+                %% this option is obsolete
+                obj.Clean=xor(obj.Clean,true); %togle state 
+            elseif ~islogical(varargin{1})
+                id = 'BSM:SetClean:InvalidType';
+                msg= 'Value passed to obj.Clean is not of type logical';
+                throw(MException(id,msg));
+            else
+                obj.Clean=varargin{1};
+            end           
+        end
+        function set.Prefix(obj,val)
+            if ischar(val)
+                obj.Prefix=val;
+            else
+                id='BSM:SetPrefix:InvalidInput';
+                msg='Input needs to be a char array';
+                throw(MException(id,msg));
+            end
+        end
+        function set.Type(obj,val)
+            %check allowed input
+            if ~isnumeric(val)|| any(~ismember(val,1:5))
+                id='BSM:Type:invalid';
+                msg='use a value in the range 1..5';
+                throw(MException(id,msg));
+            end
+            obj.Type=val;
+        end
+        function set.Acfl(obj,val)
+            %allow clearing
+            if isempty(val)
+                obj.Acfl=[];
+            end
+            if ~isnumeric(val) || any(val<1) || any(~isfinite(val))
+                id='BSM:Acfl:invalid';
+                msg='set an array of positive integer values to the Acfl field';
+                throw(MException(id,msg));
+            end
+            obj.Acfl=val;
+            
+        end
+        function set.SPMs(obj,val)
+            %set the SPMs to be analysed
+            %obj.SPMs=<val>
+            %possible values for <val>
+            %   {}      : clear the SPMs field
+            %   {'ui'}  : ask user input
+            if ~iscell(val)
+                fprintf('%s\n',...
+                    'Use any of the following options',...
+                    'obj.SPMs={}',...
+                    'obj.SPMs={''ui''}',...
+                    'obj.SPMs={<Ful path to SPM.mat>}');
+                return
+            end
+            if isempty(val) %allow empty input
+                return;
+            end
+            if strcmpi(val{1},'ui')
+                [P,sts]=spm_select([1 Inf],'mat','select the SPM.mat files',obj.SPMs,pwd,'^SPM.mat');
+                if ~logical(sts);return;end
+                if iscell(P)
+                    val=P;
+                else
+                    val=cellstr(P);
+                end
+            end
+            % capture file not found somewhat clever
+            TF=cellfun(@exist,val);
+            if any(not(TF))
+                id='BSM:SetSPMs:FileNotFound';
+                msg=sprintf('file not found\t%s\n',val{not(TF)});
+                throw(MException(id,msg));
+            end
+            obj.SPMs=val;
+        end
+    end
+    methods (Access = protected)%private functions
+        function CoreCalc(obj,idx)
+            fprintf('working on :%d: %s\n',idx,obj.SPMs{idx});
+            % how to store the current status of warning and restore curent
+            WarnStat=warning;
+            warning('OFF', 'MATLAB:MKDIR:DirectoryExists');
+            %this try catch block will handle any errors that might ocure
+            %warning will not remain off by accident
+            try
+                SpmFile=obj.SPMs{idx};
+                load(SpmFile); %#ok<LOAD>
+                %% exception block
+                %NB these exception has not been tested properly by the
+                if (SPM.xBF.order ~=1 ) %#ok<*NODEF>
+                    %TestBSMClass program
+                    %throw an error. The calling function should properly
+                    %handle the error.
+                    id='BSM:CoreCalc:SPMxBFOrder:invalid';
+                    msg=sprintf('%s\n','Problem found', ...
+                        'No canonical HRF was used to model the data', ...
+                        'Current file will be skipped');
+                    throw(MException(id,msg));
+                end
+                if (length(SPM.Sess) > 1 )
+                    id='BSM:CoreCalc:SPMSess:invalid';
+                    msg=sprintf('%s\n',...
+                        'Problem found', ...
+                        'Multiple sessions in one SPM.mat not supported', ...
+                        'Current file will be skipped');
+                    throw(MException(id,msg));
+                end
+                % end exception block
+                %% set local variables
+                SPMOrig = SPM;
+                UOrig = SPM.Sess(1).U;
+                [ActCol,~,~]=obj.GetColDef(UOrig);
+                LocalTargetFolder = fullfile(fileparts(SpmFile),'BSM');%RR replace BSM by a obj.BSMDir so user can modify this
+                %% create directory
+                mkdir(LocalTargetFolder);
+                % loop over ActCol
+                for c=1:length(ActCol)
+                    % get current col, get other col
+                    CurCol=ActCol(c);
+                    OthCol=setdiff(1:length(ActCol),CurCol);
+                    TargetFolderCurCol=fullfile(LocalTargetFolder,sprintf('col_%04d', CurCol));
+                    NEvents=numel(UOrig(CurCol).ons);
+                    DirTrialList=cell(NEvents,1);
+                    BetaTrialList=cell(NEvents,1);
+                    for e=1:NEvents%      loop over events
+                    	U=obj.GetU(UOrig,e,CurCol,OthCol);%build U
+                    	DirTrialList{e}=obj.GenSPM(U,SPMOrig,TargetFolderCurCol,e);%build SPM
+                       % obj.RunJob(DirTrialList{e});%process SPM
+                        %BetaTrialList{e}=spm_select('ExtFPList',DirTrialList{e},'^beta_0001.*.nii$',1);
+                    end%      end loop
+                    try 
+                        obj.PackBetas(BetaTrialList,LocalTargetFolder,CurCol)%pack per event result
+                        obj.CleanUp(TargetFolderCurCol);
+                    catch ME
+                        keyboard; %for debuging purpose only
+                    end
+                end % end loop act. col
+                % restore warning status
+                warning(WarnStat);
+            catch ME %turn matlab warning back on
+                warning(WarnStat);
+                rethrow(ME);
+            end
+        end
+        function CleanUp(obj,D)
+            if obj.Clean
+                try
+                    rmdir(D,'s');
+                catch ME
+                    disp(ME)
+                    id = 'BSM:Clean:Failed';
+                    msg= 'Cleaning failed';
+                    throw(MException(id,msg));
+                end
+            end
+        end
+        function PackBetas(obj,P,wd,Col)
+            FPFout=fullfile(wd,sprintf('%s_%04d.nii',obj.Prefix,Col));
+            matlabbatch{1}.spm.util.cat.vols = P;
+            matlabbatch{1}.spm.util.cat.name = FPFout;
+            matlabbatch{1}.spm.util.cat.dtype = 4;
+            spm_jobman('run', matlabbatch);
+        end
+
+        function RunJob(~,swd)
+            %% estimate model:
+            matlabbatch{1}.spm.stats.fmri_est.spmmat = {fullfile(swd,'SPM.mat')};
+            matlabbatch{1}.spm.stats.fmri_est.method.Classical = 1;
+            spm_jobman('run', matlabbatch);
+        end
+        function swd=GenSPM(obj,U,SPMOrig,TargetFolderCol,Event)
+            SPM=SPMOrig;
+            TargetFolderColTrial=fullfile(TargetFolderCol,sprintf('trial_%04d',Event));
+            mkdir(TargetFolderColTrial);
+            SPM.swd = TargetFolderColTrial;
+            SPM.Sess.U=U;
+            U=spm_get_ons(SPM,1); % to have SPM autocomplete the U structure
+            [X,~,Fc]=spm_Volterra(U,SPM.xBF.bf,SPM.xBF.Volterra);
+            
+            k=SPM.nscan(1);
+            fMRI_T  = SPM.xBF.T;
+            fMRI_T0 = SPM.xBF.T0;
+            
+            X = X((0:(k-1))*fMRI_T +fMRI_T0 + 32,:);%@JB what does the +32 do?
+            %orthogonalize
+            for k=1:numel(Fc)
+                X(:,Fc(k).i)=spm_orth(X(:,Fc(k).i));
+            end
+            
+            Xname=cell(1,size(X,2));
+            for k = 1:size(X,2)
+                Xname{k} = [sprintf('Sn(%i) ',1) SPM.Sess(1).U(k).name{1}];
+            end
+
+            %% reinstate motion parameters
+            if (obj.KeepMP)&& ~isempty(SPMOrig.Sess(1).C.C)
+                C = SPMOrig.Sess(1).C.C;
+                Cname = job.SPM_orig.Sess(1).C.name;
+                
+                X       = cat(2,X,spm_detrend(C));
+                Xname   = cat(2,Xname, Cname(:));
+            end
+            %% reinstate baseline column
+            X(:,end+1) = 1;
+            Xname = cat(2,Xname,{sprintf('Sn(%i) %s',1,'constant')});
+            %% monitor
+            if obj.MonitorState
+                figure(100);
+                clf;
+                imagesc(X);
+                pause(0.5);
+            end
+            %% Rebuild SPM
+            SPM.xX.X    = X;
+            SPM.xX.iB   = size(X,2);
+            SPM.xX.iC   = 1:size(X,2)-1;
+            SPM.xX.name = Xname;
+            try
+                SPM.xX = rmfield(SPM.xX,'W');%remove prewhitening matrix
+            catch
+            end
+            save(fullfile(SPM.swd,'SPM.mat'),'SPM');
+            swd=SPM.swd;
+        end
+        function U=GetU(obj,UOrig,Event,Col,OCol)
+            U=struct('name',{},'ons', [], 'dur', [],'orth', [], 'P', [], 'dt', [], 'u', [],'pst',[]);
+            idx=1;
+            %split out current event as first condition
+            U(idx).name{1}  = sprintf('%s_%d',UOrig(Col).name{1},Event);
+            U(idx).ons      = UOrig(Col).ons(Event);
+            U(idx).dur      = UOrig(Col).dur(Event);
+            U(idx).orth     = UOrig(Col).orth;
+            U(idx).P        = UOrig(Col).P;
+            U(idx).dt       = UOrig(Col).dt;
+            % u
+            % pst
+            %optional add all other events as the next condition
+            switch obj.Type
+                case {1,2} 
+                    idx = idx +1;
+                    OtherTrials=setdiff(1:numel(UOrig(Col).ons),Event);
+                    U(idx).name{1}  = sprintf('%s_oe',UOrig(Col).name{1});
+                    U(idx).ons      = UOrig(Col).ons(OtherTrials);
+                    U(idx).dur      = UOrig(Col).dur(OtherTrials);
+                    U(idx).orth     = UOrig(Col).orth;
+                    U(idx).P        = UOrig(Col).P;
+                    U(idx).dt       = UOrig(Col).dt;
+                    % u
+            % pst
+                case {3,4}
+                otherwise
+            end
+            % deal with other conditions
+            if ~isempty(OCol) 
+                switch obj.Type
+                    case {1,4} %put all other conditions in one column
+                        idx = idx + 1;
+                        U(idx).name{1}  = 'AllOtherConditions';
+                        U(idx).ons      = cat(1,UOrig(OCol).ons);
+                        for k = OCol
+                            if numel(UOrig(k).dur)==1
+                                dur=repmat(UOrig(k).dur,size(UOrig(k).ons));
+                            elseif numel(UOrig(k).dur)==numel(UOrig(k).ons)
+                                dur=UOrig(k).dur;
+                            else
+                                id='BSM:GetU:InvalidNumberDuration';
+                                msg=sprintf('%s\n'...
+                                    ,'Number of durations is neither one nor the number of onsets'...
+                                    );
+                                throw(MException(id,msg));
+                            end
+                            U(idx).dur=cat(1,U(idx).dur,dur);
+                            U(idx).orth     = UOrig(k).orth;
+                            U(idx).P        = UOrig(k).P;
+                            U(idx).dt       = UOrig(k).dt;
+                        end
+                    case {2,3} %keep other conditions each in their own column
+                        U=cat(2,U,UOrig(OCol));
+                    otherwise
+                end
+                %clear out fields u and pst They will be rebuild when SPM
+                %is created
+                [U(:).pst]=deal([]);
+                [U(:).u]=deal([]);
+            end
+        end
+        function [ActCol,PasCol,AllCol]=GetColDef(obj,U)
+            AllCol=1:length(U);
+            if ~isempty(obj.Acfl)
+                PasCol = AllCol(obj.Acfl); %cf line 53 in spm_bsm.m
+            else
+                PasCol = [];
+            end
+            ActCol = setdiff(AllCol,PasCol);
+        end
+        function[val,TF]=MyGetVal(~,selector,List)
+            TF=strcmpi(selector,List);
+            if any(TF) %spmpath is set as option
+                k=find(TF, 1 );%find the first instance
+                if (k+1)>numel(List)
+                    id='BSM:getval:ReadBeyondEndInput';
+                    msg=sprintf('use name value pairs\nlast entry is a name');
+                    throw(MException(id,msg));
+                end
+                val=List{k+1};
+                TF=true;
+            else
+                val=false;
+                TF=false;
+            end
+        end
+    end
+end
+

resources/BSM_Class/TestBSMClass.m

@@ -0,0 +1,123 @@
+%% BSM code validity check.  
+% This code will check the BSM class code in a very basic way. 
+% 
+%% create an instance
+B=BSM;
+%% test initialization
+clear B
+try
+    B=BSM('type',3,'Acfl',[],'SPMs'); %#ok<*NASGU> %will create an error
+catch ME
+    if strcmpi(ME.identifier,'BSM:getval:ReadBeyondEndInput')
+        disp('error properly handled')
+        disp(ME);
+    else
+        rethrow(ME);
+    end
+end
+%% initialize properly
+clear B
+B=BSM('type',1,'Acfl',[],'SPMs',{});
+%% set SPMs validity check
+B.SPMs=''; %should produce instructions on how to set SPMs
+B.SPMs={'ui'}; %get input from the user
+tmp=B.SPMs; %store for future use
+%B.SPMs={'ui'}; %calling it again will allow me to edit the list
+B.SPMs={}; %clear out SPMs
+try 
+    B.SPMs={'SPM.mat';'onzin.mat'}; %provides an error if file does not exist
+catch ME
+    if strcmpi(ME.identifier,'BSM:SetSPMs:FileNotFound')
+    disp('error properly handled')
+        disp(ME);
+    else
+        rethrow(ME)
+    end 
+end
+%%
+try 
+    B.SPMs={'SPM.mat';'onzin.mat'}; %provides an error if file does not exist
+catch ME
+    if strcmpi(ME.identifier,'BSM:SetSPMs:FileNotFound')
+    disp('error properly handled')
+        disp(ME);
+    else
+        rethrow(ME)
+    end 
+end
+B.SPMs=tmp;
+%% set Prefix check
+try
+    B.Prefix=""; %gives an error
+catch ME
+    if strcmpi(ME.identifier,'BSM:SetPrefix:InvalidInput')
+        disp('error properly handled')
+        disp(ME);
+    else
+        rethrow(ME)
+    end
+end
+B.Prefix='bs'
+%% set Clean check
+try
+    B.Clean=""; %gives an error
+catch ME
+    if strcmpi(ME.identifier,'BSM:SetClean:InvalidType')
+        disp('error properly handled')
+        disp(ME);
+    else
+        rethrow(ME)
+    end
+end
+disp(B.Clean)
+B.Clean=true;
+%% Check Type set
+try 
+    B.Type=0; %gives an error
+catch ME
+     if strcmpi(ME.identifier,'BSM:Type:invalid')
+        disp('error properly handled')
+        disp(ME);
+    else
+        rethrow(ME)
+     end
+end
+%% 
+B.Type=1; % correct call
+%% set Acfl check
+try
+    B.Acfl=0;
+catch ME
+    if strcmpi(ME.identifier,'BSM:Acfl:invalid')
+        disp('error properly handled')
+        disp(ME);
+    else
+        rethrow(ME)
+    end
+end
+%% 
+try
+    B.Acfl=-5;
+catch ME
+    if strcmpi(ME.identifier,'BSM:Acfl:invalid')
+        disp('error properly handled')
+        disp(ME);
+    else
+        rethrow(ME)
+    end
+end
+%%
+try
+    B.Acfl='test';
+catch ME
+    if strcmpi(ME.identifier,'BSM:Acfl:invalid')
+        disp('error properly handled')
+        disp(ME);
+    else
+        rethrow(ME)
+    end
+end
+%%
+B.Acfl=[];
+%% Test Calc routine
+B.Calc

resources/project/Project.xml

@@ -1,2 +0,0 @@
-<?xml version='1.0' encoding='UTF-8'?>
-<Info MetadataType="distributed" />

resources/project/ProjectData.type.Info.xml

@@ -1,2 +0,0 @@
-<?xml version='1.0' encoding='UTF-8'?>
-<Info Name="fMRI_Connectivity" />

resources/project/Root.type.Categories/FileClassCategory.type.Category.xml

@@ -1,2 +0,0 @@
-<?xml version='1.0' encoding='UTF-8'?>
-<Info ReadOnly="1" SingleValued="1" DataType="None" Name="Classification" />

resources/project/Root.type.Categories/FileClassCategory.type.Category/artifact.type.Label.xml

@@ -1,2 +0,0 @@
-<?xml version='1.0' encoding='UTF-8'?>
-<Info ReadOnly="READ_ONLY" Name="Artifact" />

resources/project/Root.type.Categories/FileClassCategory.type.Category/convenience.type.Label.xml

@@ -1,2 +0,0 @@
-<?xml version='1.0' encoding='UTF-8'?>
-<Info ReadOnly="READ_ONLY" Name="Convenience" />

resources/project/Root.type.Categories/FileClassCategory.type.Category/derived.type.Label.xml

@@ -1,2 +0,0 @@
-<?xml version='1.0' encoding='UTF-8'?>
-<Info ReadOnly="READ_ONLY" Name="Derived" />

resources/project/Root.type.Categories/FileClassCategory.type.Category/design.type.Label.xml

@@ -1,2 +0,0 @@
-<?xml version='1.0' encoding='UTF-8'?>
-<Info ReadOnly="READ_ONLY" Name="Design" />

resources/project/Root.type.Categories/FileClassCategory.type.Category/none.type.Label.xml

@@ -1,2 +0,0 @@
-<?xml version='1.0' encoding='UTF-8'?>
-<Info ReadOnly="READ_ONLY" Name="None" />

resources/project/Root.type.Categories/FileClassCategory.type.Category/other.type.Label.xml

@@ -1,2 +0,0 @@
-<?xml version='1.0' encoding='UTF-8'?>
-<Info ReadOnly="READ_ONLY" Name="Other" />

resources/project/Root.type.Categories/FileClassCategory.type.Category/test.type.Label.xml

@@ -1,2 +0,0 @@
-<?xml version='1.0' encoding='UTF-8'?>
-<Info ReadOnly="READ_ONLY" Name="Test" />

resources/project/Root.type.ProjectPath/8f05604a-e5f3-4d08-8e1e-e6e701ad05fe.type.Reference.xml

@@ -1,2 +0,0 @@
-<?xml version='1.0' encoding='UTF-8'?>
-<Info Ref="" Type="Relative" />

resources/project/uuid-3bb1211c-321b-4f32-bc29-d38252561dca.xml

@@ -1,2 +0,0 @@
-<?xml version='1.0' encoding='UTF-8'?>
-<Info />