From 4add96e8bbefb8ff5b9dc94fb51c0de0369c63c0 Mon Sep 17 00:00:00 2001
From: p147685 <p147685@XAP021.workspace.rug.nl>
Date: Mon, 30 Mar 2020 07:29:01 +0200
Subject: [PATCH]

---
 resources/BSM_Class/BSM.m          | 213 +++++++++++++++++++++++++----
 resources/BSM_Class/TestBSMClass.m |  41 ++++--
 2 files changed, 218 insertions(+), 36 deletions(-)

diff --git a/resources/BSM_Class/BSM.m b/resources/BSM_Class/BSM.m
index 6f5cca0..40f4aca 100644
--- a/resources/BSM_Class/BSM.m
+++ b/resources/BSM_Class/BSM.m
@@ -5,9 +5,12 @@ classdef BSM<handle
     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=1; %BSM type 1..4 5: none
+        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 %private properties
+    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
@@ -16,16 +19,34 @@ classdef BSM<handle
         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 %constants
+    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
@@ -97,6 +118,35 @@ classdef BSM<handle
         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))
@@ -155,7 +205,7 @@ classdef BSM<handle
             obj.SPMs=val;
         end
     end
-    methods %private functions
+    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
@@ -168,7 +218,7 @@ classdef BSM<handle
                 load(SpmFile); %#ok<LOAD>
                 %% exception block
                 %NB these exception has not been tested properly by the
-                if (SPM.xBF.order ~=1 )
+                if (SPM.xBF.order ~=1 ) %#ok<*NODEF>
                     %TestBSMClass program
                     %throw an error. The calling function should properly
                     %handle the error.
@@ -190,7 +240,7 @@ classdef BSM<handle
                 %% set local variables
                 SPMOrig = SPM;
                 UOrig = SPM.Sess(1).U;
-                [ActCol,PasCol,AllCol]=obj.GetColDef(UOrig);
+                [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);
@@ -200,15 +250,17 @@ classdef BSM<handle
                     CurCol=ActCol(c);
                     OthCol=setdiff(1:length(ActCol),CurCol);
                     TargetFolderCurCol=fullfile(LocalTargetFolder,sprintf('col_%04d', CurCol));
-                    for e=1:numel(UOrig(CurCol).ons)%      loop over events
-                    	U=obj.GetU(UOrig,e,c)%build U
-                    %          build SPM
-                    %          process SPM
-                    %          Cleanup
+                    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
-                    %      pack per event result
-                    %      Cleanup
-                    
+                    obj.PackBetas(BetaTrialList,LocalTargetFolder,CurCol)%pack per event result
+                    obj.CleanUp(TargetFolderCurCol);
                 end % end loop act. col
                 % restore warning status
                 warning(WarnStat);
@@ -217,15 +269,98 @@ classdef BSM<handle
                 rethrow(ME);
             end
         end
-        function U=GetU(obj,UOrig,Event,Col)
-            U=struct('name',{},'ons', [], 'dur',[], 'P',[], 'dt', [], 'u', [],'pst',[]);
+        function CleanUp(~,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} 
@@ -234,16 +369,46 @@ classdef BSM<handle
                     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
-            switch obj.Type
-                case {1,4} %put all other conditions in one column
-                case {2,3} %keep other conditions each in their own column
-                otherwise
+            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)
@@ -255,14 +420,6 @@ classdef BSM<handle
             end
             ActCol = setdiff(AllCol,PasCol);
         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[val,TF]=MyGetVal(~,selector,List)
             TF=strcmpi(selector,List);
             if any(TF) %spmpath is set as option
diff --git a/resources/BSM_Class/TestBSMClass.m b/resources/BSM_Class/TestBSMClass.m
index 71a920a..966b7f3 100644
--- a/resources/BSM_Class/TestBSMClass.m
+++ b/resources/BSM_Class/TestBSMClass.m
@@ -3,7 +3,7 @@
 % 
 %% create an instance
 B=BSM;
-%% 
+%% test initialization
 clear B
 try
     B=BSM('type',3,'Acfl',[],'SPMs'); %#ok<*NASGU> %will create an error
@@ -15,7 +15,7 @@ catch ME
         rethrow(ME);
     end
 end
-%% 
+%% initialize properly
 clear B
 B=BSM('type',1,'Acfl',[],'SPMs',{});
 %% set SPMs validity check
@@ -34,19 +34,44 @@ catch ME
         rethrow(ME)
     end 
 end
-B.SPMs=tmp;
-%% set Type check
-try
-    B.Type='onzin'; %gives an error
+%%
+try 
+    B.SPMs={'SPM.mat';'onzin.mat'}; %provides an error if file does not exist
 catch ME
-    if strcmpi(ME.identifier,'BSM:Type:invalid')
+    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