Fix Issue 484 --- pspm_sf

src/pspm_cfg/pspm_cfg_run_sf.m

@@ -78,7 +78,9 @@
     options.fresp = job.fresp;
 end
 if ~isempty(job.missing)
-  options.missing = job.missing;
+  if ischar(job.missing.missingdata)
+    options.missing = job.missing;
+  end
 end
 options.dispwin = job.dispwin;
 options.dispsmallwin = job.dispsmallwin;

src/pspm_cfg/pspm_cfg_sf.m

@@ -165,6 +165,8 @@
 filter.values = {filter_def, filter_edit};
 filter.help   = {'Specify how you want filter the SCR data.'};
 
+
+
 %% Epochs
 epochfile         = cfg_files;
 epochfile.name    = 'Epoch File';
@@ -191,7 +193,7 @@
 mrk_chan.name    = 'Marker Channel';
 mrk_chan.tag     = 'mrk_chan';
 mrk_chan.strtype = 'i';
-mrk_chan.val     = {0};
+mrk_chan.val     = {1};
 mrk_chan.num     = [1 1];
 mrk_chan.help    = {['Indicate the marker channel. By default the first marker channel is ' ...
     'assumed to contain the relevant markers.'], ['Markers are only used if you have ' ...
@@ -292,14 +294,29 @@
 fresp.help    = {'Frequency of responses to model.'};
 fresp.hidden  = true;
 
-missing         = cfg_files;
-missing.name    = 'Missing epoch file';
-missing.tag     = 'missingfile';
-missing.num     = [1 1];
-missing.filter  = '.*\.(mat|MAT)$';
-missing.help    = {['Missing (e.g. artefact) epochs in the data file, where ',...
+
+
+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.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.']};
 
+missing        = cfg_choice;
+missing.name   = 'Missing Epoch Settings';
+missing.tag    = 'missing';
+missing.val    = {missingepoch_none};
+missing.values = {missingepoch_none, missingepoch_file};
+missing.help   = {'Specify whether you would like to include missing epochs.'};
+
 % Show figures
 dispwin         = cfg_menu;
 dispwin.name    = 'Display Progress Window';

src/pspm_dcm_inv.m

@@ -553,9 +553,10 @@
       priors.SigmaX0 = zeros(7);
       priors.muX0 = zeros(7, 1);
       if trl == 1
-        priors.muX0(1) = mean(y(1:3));%) - min(y);
-        priors.muX0(2) = mean(diff(y(1:3)));
-        priors.muX0(3) = diff(diff(y(1:3)));
+        y_non_nan = y(~isnan(y));
+        priors.muX0(1) = mean(y_non_nan(1:3));%) - min(y);
+        priors.muX0(2) = mean(diff(y_non_nan(1:3)));
+        priors.muX0(3) = diff(diff(y_non_nan(1:3)));
         priors.muX0(7) = 0;%min(y);
         for n = [1:3 7]
           priors.SigmaX0(n, n) = 1e-2;

src/pspm_options.m

@@ -411,6 +411,7 @@
     options = autofill(options,'dispwin',                 1,          0                 );
     options = autofill(options,'fresp',                   0.5,        '>=', 0           );
     options = autofill(options,'marker_chan_num',         1,          '*Int*Char'       );
+    options = autofill(options,'missingthresh',           2,          '>', 0            );
     options = autofill(options,'overwrite',               1,          0                 );
     options = autofill(options,'threshold',               0.1,        '>', 0            );
     options = autofill(options,'missingthresh',           2,          '>', 0            );
@@ -579,7 +580,7 @@
               if contains(optional_value, '*Int')
                 flag_is_allowed_value = flag_is_allowed_value || ...
                   all([isnumeric(options.(field_name)), ...
-                  options.(field_name)>0, ...
+                  options.(field_name)>=0, ...
                   mod(options.(field_name), 1)==0]);
               end
               if contains(optional_value, '*Struct')

src/pspm_overwrite.m

@@ -96,5 +96,7 @@
   case 2
     varargout{1} = sts;
     varargout{2} = overwrite_final;
+  otherwise
+    varargout{1} = overwrite_final;
 end
 return

src/pspm_sf.m

@@ -35,7 +35,7 @@
 %   ┌─────────options
 %   ├──────.overwrite:  [logical] [default: determined by pspm_overwrite]
 %   │                   Define whether to overwrite existing output files or not.
-%   ├.marker_chan_num:  [integer] 
+%   ├.marker_chan_num:  [integer]
 %   │                   marker channel number
 %   │                   if undefined or 0, first marker channel is used.
 %   │ * Additional options for individual methods:
@@ -77,7 +77,7 @@
 outfile = [];
 sts = -1;
 %% 2 Check input
-% 2.1 Check missing input
+% 2.1 Check missing input --
 if nargin<1
   warning('ID:invalid_input', 'Nothing to do.'); return;
 elseif nargin<2
@@ -92,7 +92,7 @@
 elseif ~isfield(model, 'timing') && ~strcmpi(model.timeunits, 'file')
   warning('ID:invalid_input', 'No epochs specified.'); return;
 end
-% 2.2 Check faulty input
+% 2.2 Check faulty input --
 if ~ischar(model.datafile) && ~iscell(model.datafile)
   warning('ID:invalid_input', 'Input data must be a cell or string.'); return;
 elseif ~ischar(model.modelfile) && ~iscell(model.modelfile)
@@ -115,15 +115,15 @@
 if ischar(model.modelfile)
   model.modelfile = {model.modelfile};
 end
-% 2.4 Check number of files
+% 2.4 Check number of files --
 if ~strcmpi(model.timeunits, 'whole') && numel(model.datafile) ~= numel(model.timing)
   warning('ID:number_of_elements_dont_match',...
     'Number of data files and epoch definitions does not match.'); return;
 elseif numel(model.datafile) ~= numel(model.modelfile)
   warning('ID:number_of_elements_dont_match',...
     'Number of data files and model files does not match.'); return;
 end
-% 2.5 check methods
+% 2.5 check methods --
 if ~isfield(model, 'method')
   model.method = {'dcm'};
 elseif ischar(model.method)
@@ -162,31 +162,31 @@
     end
   end
 end
-% 2.6 Check timing
+% 2.6 Check timing --
 if strcmpi(model.timeunits, 'whole')
   epochs = repmat({[1 1]}, numel(model.datafile), 1);
 else
-  for iSn = 1:numel(model.datafile)
-    [sts_get_timing, epochs{iSn}] = pspm_get_timing('epochs', model.timing{iSn}, model.timeunits);
+  for iFile = 1:numel(model.datafile)
+    [sts_get_timing, epochs{iFile}] = pspm_get_timing('epochs', model.timing{iFile}, model.timeunits);
     if sts_get_timing == -1
       warning('ID:invalid_input', 'Call of pspm_get_timing failed.');
       return;
     end
   end
 end
-% 2.7 Check filter
+% 2.7 Check filter --
 if ~isfield(model, 'filter')
   model.filter = settings.dcm{2}.filter;
 elseif ~isfield(model.filter, 'down') || ~isnumeric(model.filter.down)
   warning('ID:invalid_input', 'Filter structure needs a numeric ''down'' field.'); return;
 end
-% 2.8 Set options
+% 2.8 Set options --
 try model.channel; catch, model.channel = 'scr'; end
 options = pspm_options(options, 'sf');
 if options.invalid
   return
 end
-% 2.9 Set missing epochs
+% 2.9 Set missing epochs --
 if ~isfield(model, 'missing')
   model.missing = cell(numel(model.datafile), 1);
 elseif ischar(model.missing) || isnumeric(model.missing)
@@ -197,14 +197,15 @@
   return
 end
 %% 3 Get data
+missing = cell(size(model.missing));
 for iFile = 1:numel(model.datafile)
-  % 3.1 User output
+  % 3.1 User output --
   fprintf('SF analysis: %s ...', model.datafile{iFile});
-  % 3.2 Check whether model file exists
+  % 3.2 Check whether model file exists --
   if ~pspm_overwrite(model.modelfile, options)
     return
   end
-  % 3.3 get and filter data
+  % 3.3 get and filter data --
   [sts_load_data, ~, data] = pspm_load_data(model.datafile{iFile}, model.channel);
   if sts_load_data < 0, return; end
   Y{1} = data{1}.data; sr(1) = data{1}.header.sr;
@@ -214,23 +215,23 @@
     warning('ID:invalid_input', 'Call of pspm_prepdata failed.');
     return;
   end
-  % 3.4 Check data units
+  % 3.4 Check data units --
   if ~strcmpi(data{1}.header.units, 'uS') && any(strcmpi('dcm', method))
     fprintf(['\nYour data units are stored as %s, ',...
       'and the method will apply an amplitude threshold in uS. ',...
       'Please check your results.\n'], ...
       data{1}.header.units);
   end
-  % 3.5 Get missing epochs
+  % 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');
-    % sort missing epochs
+    [~, 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{iSn}, 1)
+      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, :) = [];
@@ -240,33 +241,16 @@
   else
     missing{iFile} = [];
   end
-  % 3.6 Get marker data
+  % 3.6 Get marker data --
   if any(strcmp(model.timeunits, {'marker', 'markers'}))
     if options.marker_chan_num
-      [nsts, ~, ndata] = pspm_load_data(model.datafile, options.marker_chan_num);
-      if nsts == -1
-        warning('ID:invalid_input', 'Could not load data');
-        return;
-      end
-      if ~strcmp(ndata{1}.header.chantype, 'marker')
-        warning('ID:invalid_option', ...
-          ['Channel %i is no marker channel. ',...
-          'The first marker channel in the file is used instead'],...
-          options.marker_chan_num);
-        [nsts, ~, ~] = pspm_load_data(model.datafile, 'marker');
-        if nsts == -1
-          warning('ID:invalid_input', 'Could not load data');
-          return;
-        end
-      end
+      [nsts, ~, ndata] = pspm_load_data(model.datafile{iFile}, options.marker_chan_num);
+      if nsts == -1; warning('ID:invalid_input', 'Could not load data'); return; end
     else
-      [nsts, ~, ~] = pspm_load_data(model.datafile, 'marker');
-      if nsts == -1
-        warning('ID:invalid_input', 'Could not load data');
-        return;
-      end
+      [nsts, ~, ndata] = pspm_load_data(model.datafile{iFile}, 'marker');
+      if nsts == -1; warning('ID:invalid_input', 'Could not load data'); return; end
     end
-    events = data{1}.data;
+    events = ndata{1}.data;
   end
   for iEpoch = 1:size(epochs{iFile}, 1)
     if iEpoch > 1, fprintf('\n\t\t\t'); end
@@ -279,7 +263,7 @@
         case 'samples'
           win = round(epochs{iFile}(iEpoch, :) * sr(datatype(k)) / sr(1));
         case 'markers'
-          win = round(events(epochs{1}(iEpoch, :)) * sr(datatype(k)));
+          win = round(events(epochs{iFile}(iEpoch, :)) * sr(datatype(k)));
         case 'whole'
           win = [1 numel(Y{datatype(k)})];
       end
@@ -290,7 +274,7 @@
         win(1) = max(win(1), 1);
         win(2) = min(win(2), numel(Y{datatype(k)}));
       end
-      % 3.6.1 collect information
+      % 3.6.1 collect information --
       sf.model{k}(iEpoch).modeltype = method{k};
       sf.model{k}(iEpoch).boundaries = squeeze(epochs{iFile}(iEpoch, :));
       sf.model{k}(iEpoch).timeunits  = model.timeunits;
@@ -299,8 +283,17 @@
       %
       escr = Y{datatype(k)}(win(1):win(end));
       sf.model{k}(iEpoch).data = escr;
-      % 3.6.2 do the analysis and collect results
-      model_analysis = struct('scr', escr, 'sr', sr(datatype(k)));
+      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);
+      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;

src/pspm_sf_auc.m

@@ -1,4 +1,4 @@
-function varargout = pspm_sf_auc(scr, sr, options)
+function varargout = pspm_sf_auc(model, options)
 % ● Description
 %   pspm_sf_auc returns the integral/area under the curve of an SCR time series
 % ● Format
@@ -28,6 +28,10 @@
 if nargin < 1
   warning('No data specified'); return;
 end;
+try model.scr; catch, warning('Input data is not defined.'); return; end
+try model.sr; catch, warning('Sample rate is not defined.'); return; end
+scr = model.scr;
+sr = model.sr;
 scr = scr - min(scr);
 auc = mean(scr);
 sts = 1;