Fix Issue 484 --- pspm_sf

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,'theta',                   [0.923581, ...

src/pspm_sf.m

@@ -71,7 +71,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
@@ -86,7 +86,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)
@@ -109,15 +109,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)
@@ -156,7 +156,7 @@
     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
@@ -168,19 +168,19 @@
     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)
@@ -191,15 +191,16 @@
   return
 end
 %% 3 Get data
-for iFile = 1:numel(model.datafile)
-  % 3.1 User output
-  fprintf('SF analysis: %s ...', model.datafile{iFile});
-  % 3.2 Check whether model file exists
+missing = cell(size(model.missing));
+for iSn = 1:numel(model.datafile)
+  % 3.1 User output --
+  fprintf('SF analysis: %s ...', model.datafile{iSn});
+  % 3.2 Check whether model file exists --
   if ~pspm_overwrite(model.modelfile, options)
     return
   end
-  % 3.3 get and filter data
-  [sts_load_data, ~, data] = pspm_load_data(model.datafile{iFile}, model.channel);
+  % 3.3 get and filter data --
+  [sts_load_data, ~, data] = pspm_load_data(model.datafile{iSn}, model.channel);
   if sts_load_data < 0, return; end
   Y{1} = data{1}.data; sr(1) = data{1}.header.sr;
   model.filter.sr = sr(1);
@@ -208,33 +209,33 @@
     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
-  if ~isempty(model.missing{iFile})
-    [~, missing{iFile}] = pspm_get_timing('epochs', ...
-      model.missing{iFile}, 'seconds');
-    % sort missing epochs
-    if size(missing{iFile}, 1) > 0
-      [~, sortindx] = sort(missing{iFile}(:, 1));
-      missing{iFile} = missing{iFile}(sortindx,:);
+  % 3.5 Get missing epochs --
+  % 3.5.1 Load missing epochs --
+  if ~isempty(model.missing{iSn})
+    [~, missing{iSn}] = pspm_get_timing('epochs', model.missing{iSn}, 'seconds');
+  % 3.5.2 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{iFile}(k, 1) <= missing{iFile}(k - 1, 2)
-          missing{iFile}(k, 1) =  missing{iFile}(k - 1, 1);
-          missing{iFile}(k - 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{iFile} = [];
+    missing{iSn} = [];
   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);
@@ -262,39 +263,47 @@
     end
     events = data{1}.data;
   end
-  for iEpoch = 1:size(epochs{iFile}, 1)
+  for iEpoch = 1:size(epochs{iSn}, 1)
     if iEpoch > 1, fprintf('\n\t\t\t'); end
     fprintf('epoch %01.0f ...', iEpoch);
     for k = 1:numel(method)
       fprintf('%s ', method{k});
       switch model.timeunits
         case 'seconds'
-          win = round(epochs{iFile}(iEpoch, :) * sr(datatype(k)));
+          win = round(epochs{iSn}(iEpoch, :) * sr(datatype(k)));
         case 'samples'
-          win = round(epochs{iFile}(iEpoch, :) * sr(datatype(k)) / sr(1));
+          win = round(epochs{iSn}(iEpoch, :) * sr(datatype(k)) / sr(1));
         case 'markers'
           win = round(events(epochs{1}(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});
+        warning('\nEpoch %2.0f outside of file %s ...', iEpoch, model.modelfile{iSn});
       else
         % correct issues with using 'round'
         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).boundaries = squeeze(epochs{iSn}(iEpoch, :));
       sf.model{k}(iEpoch).timeunits  = model.timeunits;
       sf.model{k}(iEpoch).samples    = win;
       sf.model{k}(iEpoch).sr         = sr(datatype(k));
       %
       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{iSn})
+        model.missing_data = zeros(size(escr));
+        model.missing_data((missing{iSn}(:,1)+1):(missing{iSn}(:,2)+1)) = 1;
+      end
+      % 3.6.2 do the analysis and collect results --
+      if any(missing{iSn})
+        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;
@@ -308,16 +317,16 @@
   end
   sf.names = method(:);
   sf.infos.date = date;
-  sf.infos.file = model.modelfile{iFile};
-  sf.modelfile = model.modelfile{iFile};
+  sf.infos.file = model.modelfile{iSn};
+  sf.modelfile = model.modelfile{iSn};
   sf.data = Y;
   if exist('events','var'), sf.events = events; end
   sf.input = model;
   sf.options = options;
   sf.modeltype = 'sf';
   sf.modality = settings.modalities.sf;
-  save(model.modelfile{iFile}, 'sf');
-  outfile = model.modelfile(iFile);
+  save(model.modelfile{iSn}, 'sf');
+  outfile = model.modelfile(iSn);
   fprintf('\n');
 end
 sts = 1;

src/pspm_sf_dcm.m

@@ -5,18 +5,6 @@
 %   the input data is assumed to be in mcS, and sampling rate in Hz
 % ● Format
 %   function out = pspm_sf_dcm(model, options)
-% ● Output
-%         out:  output
-%          .n:  number of responses above threshold
-%          .f:  frequency of responses above threshold in Hz
-%         .ma:  mean amplitude of responses above threshold
-%          .t:  timing of (all) responses
-%          .a:  amplitude of (all) responses
-%      .theta:  parameters used for f_SF
-%  .threshold:  threshold
-%         .if:  initial frequency for f_SF
-%       .yhat:  fitted time series
-%      .model:  information about the DCM inversion
 % ● Arguments
 %   ┌──────model
 %   │ ▶︎ Mandatory
@@ -25,7 +13,6 @@
 %   ├────────.sr:  sampling rate in Hz
 %   │ ▶︎ Optional
 %   └.missing_data: the datafile of missing epochs
-%
 %   ┌────options:  options structure
 %   ├─.threshold:  threshold for SN detection (default 0.1 mcS)
 %   ├─────.theta:  a (1 x 5) vector of theta values for f_SF
@@ -36,6 +23,18 @@
 %   │              display intermediate windows (default 0);
 %   └─.missingthresh:
 %               threshold value for controlling missing epochs (default 2s).
+% ● Output
+%   ┌────────out:  output
+%   ├─────────.n:  number of responses above threshold
+%   ├─────────.f:  frequency of responses above threshold in Hz
+%   ├────────.ma:  mean amplitude of responses above threshold
+%   ├─────────.t:  timing of (all) responses
+%   ├─────────.a:  amplitude of (all) responses
+%   ├─────.theta:  parameters used for f_SF
+%   ├─.threshold:  threshold
+%   ├────────.if:  initial frequency for f_SF
+%   ├──────.yhat:  fitted time series
+%   └─────.model:  information about the DCM inversion
 % ● References
 %   Bach DR, Daunizeau J, Kuelzow N, Friston KJ, & Dolan RJ (2011). Dynamic
 %   causal modelling of spontaneous fluctuations in skin conductance.
@@ -51,12 +50,10 @@
 end
 sts = -1;
 tstart = tic;
-
 %% 2 Check input arguments
 % 2.1 set model ---
 try model.scr; catch, warning('Input data is not defined.'); return; end
 try model.sr; catch, warning('Sample rate is not defined.'); return; end
-try model.missing_data; catch, warning('Missing data file is not defined.'); end
 % 2.2 Validate parameters ---
 if ~isnumeric(model.sr) || numel(model.sr) > 1
   errmsg = sprintf('No valid sample rate given.');
@@ -115,7 +112,8 @@
 y = model.scr;
 y = y - min(y);
 % 4.4 determine initial conditions
-x0 = y(1:3);
+y_non_nan = y(~isnan(y));
+x0 = y_non_nan(1:3);
 X0(1, 1)   = mean(x0);
 X0(2, 1)   = mean(diff(x0));
 X0(3, 1)   = diff(diff(x0));
@@ -124,14 +122,14 @@
 u = [];
 u(1, :) = (1:numel(y))/model.sr;
 u(2, :) = nresp;
-priors.muTheta = theta(1:3)';
+priors.muTheta = transpose(theta(1:3));
 priors.muTheta(4:2:(2 * nresp + 3)) = 1/fresp * (0:(nresp-1));
 priors.muTheta(5:2:(2 * nresp + 4)) = -10;
 dim.n_theta = numel(priors.muTheta);    % nb of evolution parameters
 priors.SigmaTheta = zeros(dim.n_theta);
 for k = (4:2:(2 * nresp + 3)), priors.SigmaTheta(k, k) = 1e-2;end
 for k = (5:2:(2 * nresp + 4)), priors.SigmaTheta(k, k) = 1e2; end
-priors.muPhi = phi';
+priors.muPhi = transpose(phi);
 priors.SigmaPhi = zeros(dim.n_phi);
 priors.SigmaX0 = 1e-8*eye(dim.n);
 options.priors = priors;
@@ -153,16 +151,23 @@
 end
 miss_epoch = [ymissing_start(:),ymissing_end(:)];
 flag_missing_too_long = 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.'];
+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.'];
+    flag_missing_too_long = 1;
+  else
+    warning_message = ['Imported data includes miss epoches (over ',...
+      num2str(options.missingthresh), 's), but the trial has been allowed. ',...
+      'Please adjust options.missingthresh to skip if you wish.'];
+  end
   warning('ID:missing_data', warning_message);
-  flag_missing_too_long = 1;
 end
 options.isYout = ymissing(:)';
 %% 5 Extract parameters
 if ~flag_missing_too_long
-  [posterior, output] = VBA_NLStateSpaceModel(y(:)',u,f_fname,g_fname,dim,options);
+  [~,y_interpolated] = pspm_interpolate(y,struct());
+  [posterior, output] = VBA_NLStateSpaceModel(y_interpolated(:)',u,f_fname,g_fname,dim,options);
   for i = 1:length(output)
     output(i).options = rmfield(output(i).options, 'hf');
   end