Merge pull request #513 from bachlab/fix-issue-499

fix issue 499

src/pspm_cfg/pspm_cfg_run_sf.m

@@ -78,8 +78,8 @@
     options.fresp = job.fresp;
 end
 if ~isempty(job.missing)
-  if ischar(job.missing.missingdata)
-    options.missing = job.missing;
+  if ischar(job.missing.missingepoch_include.missingepoch_file{1})
+    model.missing = job.missing.missingepoch_include.missingepoch_file{1};
   end
 end
 options.dispwin = job.dispwin;

src/pspm_cfg/pspm_cfg_sf.m

@@ -296,25 +296,33 @@
 
 
 
-missingepoch_none        = cfg_const;
-missingepoch_none.name   = 'Not added';
-missingepoch_none.tag    = 'missingdata';
-missingepoch_none.val    = {0};
-missingepoch_none.help   = {'Do not add missing epochs.'};
+
 
 missingepoch_file         = cfg_files;
 missingepoch_file.name    = 'Missing epoch file';
-missingepoch_file.tag     = 'missingdata';
+missingepoch_file.tag     = 'missingepoch_file';
 missingepoch_file.num     = [1 1];
 missingepoch_file.filter  = '.*\.(mat|MAT)$';
 missingepoch_file.help    = {['Missing (e.g. artefact) epochs in the data file, where ',...
                   'data must always be specified in seconds.']};
 
+missingepoch_none        = cfg_const;
+missingepoch_none.name   = 'Do not add';
+missingepoch_none.tag    = 'missingepoch_none';
+missingepoch_none.val    = {0};
+missingepoch_none.help   = {'Do not add missing epochs.'};
+
+missingepoch_include         = cfg_branch;
+missingepoch_include.name    = 'Add';
+missingepoch_include.tag     = 'missingepoch_include';
+missingepoch_include.val     = {missingepoch_file};
+missingepoch_include.help    = {'Add missing epoch file'};
+
 missing        = cfg_choice;
 missing.name   = 'Missing Epoch Settings';
 missing.tag    = 'missing';
 missing.val    = {missingepoch_none};
-missing.values = {missingepoch_none, missingepoch_file};
+missing.values = {missingepoch_none, missingepoch_include};
 missing.help   = {'Specify whether you would like to include missing epochs.'};
 
 % Show figures

src/pspm_dcm.m

@@ -339,20 +339,8 @@
 
   % load and check existing missing data (if defined)
   if ~isempty(model.missing{iSn})
-    [~, missing{iSn}] = pspm_get_timing('epochs', ...
+    [~, missing{iSn}] = pspm_get_timing('missing', ...
       model.missing{iSn}, 'seconds');
-    % sort missing epochs
-    if size(missing{iSn}, 1) > 0
-      [~, sortindx] = sort(missing{iSn}(:, 1));
-      missing{iSn} = missing{iSn}(sortindx,:);
-      % check for overlap and merge
-      for k = 2:size(missing{iSn}, 1)
-        if missing{iSn}(k, 1) <= missing{iSn}(k - 1, 2)
-          missing{iSn}(k, 1) =  missing{iSn}(k - 1, 1);
-          missing{iSn}(k - 1, :) = [];
-        end
-      end
-    end
   else
     missing{iSn} = [];
   end

src/pspm_get_timing.m

@@ -7,6 +7,7 @@
 %   [sts, multi]  = pspm_get_timing('onsets', intiming, timeunits)
 %   [sts, epochs] = pspm_get_timing('epochs', epochs, timeunits)
 %   [sts, events] = pspm_get_timing('events', events)
+%   [sts, epochs] = pspm_get_timing('missing', epochs, timeunits)
 %   for recursive calls also:
 %   [sts, epochs] = pspm_get_timing('file', filename)
 % ● Arguments
@@ -82,7 +83,7 @@
   intiming = varargin{2};
 end
 
-if ~ismember(model, {'onsets', 'epochs', 'events', 'file'})
+if ~ismember(model, {'onsets', 'epochs', 'missing', 'events', 'file'})
   warning('ID:invalid_input', 'Invalid input. I don''t know what to do.');
 
   return;
@@ -91,7 +92,7 @@
 
 
 switch model
-  case {'onsets', 'epochs'}
+  case {'onsets', 'epochs', 'missing'}
     if nargin < 3
       warning('ID:invalid_input', 'Time units unspecified');  return;
     else
@@ -352,7 +353,7 @@
     end
 
 
-    % Epoch information for SF and GLM (model.missing)
+    % Epoch information for SF and recursive call from option "missing"
     % ------------------------------------------------------------------------
   case 'epochs'
     % get epoch information from file or from input --
@@ -380,7 +381,6 @@
       if filewarning
         warning('File %s is not a valid epochs or onsets file', ...
           intiming);
-
         return;
       end
     else
@@ -393,6 +393,10 @@
         if size(outtiming, 2) ~= 2
           warning(['Epochs must be specified by a e x 2 vector', ...
             'of onset/offsets.']);  return;
+        else
+            if any(diff(outtiming, [], 2) < 0)
+                warning('Offsets must be larger than onsets.');  return;
+            end
         end
       else
         warning('Unknown epoch definition format.');  return;
@@ -408,6 +412,26 @@
         'time units are ''%s'''], timeunits);  return;
     end
 
+       % Missing epoch information for GLM and DCM
+    % ------------------------------------------------------------------------
+  case 'missing'
+    [sts, missepochs] = pspm_get_timing('epochs', intiming, timeunits);
+     if sts < 1, return; end
+     % sort & merge missing epochs
+    if size(missepochs, 1) > 0
+      [~, sortindx] = sort(missepochs(:, 1));
+      missepochs = missepochs(sortindx,:);
+      % check for overlap and merge
+      overlapindx = zeros(size(missepochs, 1), 1);
+      for k = 2:size(missepochs, 1)
+        if missepochs(k, 1) <= missepochs(k - 1, 2)
+          missepochs(k, 1) =  missepochs(k - 1, 1);
+          overlapindx(k - 1) = 1;
+        end
+      end
+      missepochs(logical(overlapindx), :) = [];
+    end
+    outtiming = missepochs;
 
     % Event information for DCM
     % ------------------------------------------------------------------------

src/pspm_glm.m

@@ -368,7 +368,7 @@
     if isempty(model.missing{iSn})
       sts = 1; missing{iSn} = [];
     else
-      [sts, missing{iSn}] = pspm_get_timing('epochs', model.missing{iSn}, 'seconds');
+      [sts, missing{iSn}] = pspm_get_timing('missing', model.missing{iSn}, 'seconds');
     end
     if sts == -1, warning('ID:invalid_input',...
         'Failed to call pspm_get_timing'); return; end

src/pspm_interp1.m

@@ -69,11 +69,10 @@
 end
 X_body = X(index_non_nan_full(1):index_non_nan_full(end));
 % processing body
-index = 1:length(X_body);
-index_nan = index(isnan(X_body) | index_missing(index_non_nan_full(1):index_non_nan_full(end))) ;
-index_non_nan = 1 - index_nan;
-if ~isempty(index_nan)
-  X_body_interp = interp1(index_non_nan,X_body(index_non_nan),index_nan);
+index_nan = zeros(size(X_body));
+index_nan(isnan(X_body) | index_missing(index_non_nan_full(1):index_non_nan_full(end))) = 1;
+if any(index_nan)
+  X_body_interp = interp1(find(~index_nan),X_body(~index_nan), (1:numel(X_body))');
 else
   X_body_interp = X_body;
 end

src/pspm_sf.m

@@ -76,6 +76,13 @@
 end
 outfile = [];
 sts = -1;
+switch nargout
+  case 1
+    varargout{1} = outfile;
+  case 2
+    varargout{1} = sts;
+    varargout{2} = outfile;
+end
 %% 2 Check input
 % 2.1 Check missing input --
 if nargin<1
@@ -225,21 +232,11 @@
       data{end}.header.units);
   end
   % 3.5 Get missing epochs --
-  % 3.5.1 Load missing epochs --
   if ~isempty(model.missing{iFile})
-    [~, missing{iFile}] = pspm_get_timing('epochs', model.missing{iFile}, 'seconds');
-  % 3.5.2 sort missing epochs --
-    if size(missing{iFile}, 1) > 0
-      [~, sortindx] = sort(missing{iFile}(:, 1));
-      missing{iFile} = missing{iFile}(sortindx,:);
-      % check for overlap and merge
-      for k = 2:size(missing{iFile}, 1)
-        if missing{iFile}(k, 1) <= missing{iFile}(k - 1, 2)
-          missing{iFile}(k, 1) =  missing{iFile}(k - 1, 1);
-          missing{iFile}(k - 1, :) = [];
-        end
-      end
-    end
+    [~, missing{iFile}] = pspm_get_timing('missing', model.missing{iFile}, 'seconds');
+    model.missing_data = zeros(size(y{2}));
+    missing_index = pspm_time2index(missing{iFile}, sr(datatype(k)));
+    model.missing_data((missing_index(:,1)+1):(missing_index(:,2)+1)) = 1;
   else
     missing{iFile} = [];
   end
@@ -260,9 +257,9 @@
         return;
       end
     end
+    % use first marker channel
+    events{iFile} = ndata{1}.data(:);
   end
-  % use first marker channel
-  events{iFile} = data{1}.data(:);
 
 
   for iEpoch = 1:size(epochs{iFile}, 1)
@@ -276,17 +273,23 @@
         case 'samples'
           win = round(epochs{iFile}(iEpoch, :) * sr(datatype(k)) / sr(1));
         case 'markers'
-          win = round(events(epochs{iFile}(iEpoch, :)) * sr(datatype(k)));
+          win = round(events{iFile}(epochs{iFile}(iEpoch, :)) * sr(datatype(k)));
         case 'whole'
           win = [1 numel(y{datatype(k)})];
       end
       if any(win > numel(y{datatype(k)}) + 1) || any(win < 0)
         warning('\nEpoch %2.0f outside of file %s ...', iEpoch, model.modelfile{iFile});
+        inv_flag = 0;
       else
+        inv_flag = 1;
         % correct issues with using 'round'
         win(1) = max(win(1), 1);
         win(2) = min(win(2), numel(y{datatype(k)}));
       end
+      if diff(win) < 4
+          warning('\nEpoch %2.0f contains insufficient data ...', iEpoch);
+          inv_flag = 0;
+      end
       % 3.6.1 collect information --
       sf.model{k}(iEpoch).modeltype = method{k};
       sf.model{k}(iEpoch).boundaries = squeeze(epochs{iFile}(iEpoch, :));
@@ -296,23 +299,24 @@
       %
       escr = y{datatype(k)}(win(1):win(end));
       sf.model{k}(iEpoch).data = escr;
-      if any(missing{iFile})
-        model.missing_data = zeros(size(escr));
-        missing_index = pspm_time2index(missing, sr(datatype(k)));
-        model.missing_data((missing_index{iFile}(:,1)+1):(missing_index{iFile}(:,2)+1)) = 1;
-      end
+
       % 3.6.2 do the analysis and collect results --
-      if any(missing{iFile})
-        model_analysis = struct('scr', escr, 'sr', sr(datatype(k)), 'missing_data', model.missing_data);
+      if ~isempty(model.missing{iFile})
+        model_analysis = struct('scr', escr, 'sr', sr(datatype(k)), 'missing_data', model.missing_data(win(1):win(end)));
       else
         model_analysis = struct('scr', escr, 'sr', sr(datatype(k)));
       end
-      invrs = fhandle{k}(model_analysis, options);
-      if any(strcmpi(method{k}, {'dcm', 'mp'}))
-        sf.model{k}(iEpoch).inv     = invrs;
+      if inv_flag ~= 0
+        invrs = fhandle{k}(model_analysis, options);
+        sf.model{k}(iEpoch).inv = invrs;
+      else
+        sf.model{k}(iEpoch).inv = [];
+      end
+      if inv_flag == 0
+        sf.stats(iEpoch, k) = NaN;
+      elseif any(strcmpi(method{k}, {'dcm', 'mp'}))
         sf.stats(iEpoch, k)         = invrs.f;
       else
-        sf.model{k}(iEpoch).stats   = invrs;
         sf.stats(iEpoch, k)         = invrs;
       end
     end

src/pspm_sf_dcm.m

@@ -150,12 +150,13 @@
 flag_missing_too_long = 0;
 if any(diff(miss_epoch, 1, 2)/model.sr > 0)
   if any(diff(miss_epoch, 1, 2)/model.sr > options.missingthresh)
-    warning_message = ['Imported data includes too long miss epoches (over ',...
-      num2str(options.missingthresh), 's), thus estimation has been skipped.'];
+    warning_message = ['Epoch includes missing data of more than ',...
+      num2str(options.missingthresh), ' s, thus estimation has been skipped. ', ...
+      'Please adjust options.missingthresh to proceed if you wish.'];
     flag_missing_too_long = 1;
   else
-    warning_message = ['Imported data includes miss epoches (over ',...
-      num2str(options.missingthresh), 's), but the trial has been allowed. ',...
+    warning_message = ['Epoch includes missing data of less than ',...
+      num2str(options.missingthresh), ' s, hence estimation is proceeding. ', ...
       'Please adjust options.missingthresh to skip if you wish.'];
   end
   warning('ID:missing_data', warning_message);