From 868ee9dfc343eb88560b26f1bc95e2555cd771d4 Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Mon, 15 May 2023 15:57:45 +0100
Subject: [PATCH 01/48] update missthreshold for sf

---
 src/pspm_options.m |  1 +
 src/pspm_sf_dcm.m  | 14 ++++++++++----
 2 files changed, 11 insertions(+), 4 deletions(-)

diff --git a/src/pspm_options.m b/src/pspm_options.m
index 1a0daa9e1..eabde9e62 100644
--- a/src/pspm_options.m
+++ b/src/pspm_options.m
@@ -412,6 +412,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, 2]            );
     options = autofill(options,'threshold',               0.1,        '>', 0            );
     options = autofill(options,'theta',                   [0.923581, ...
diff --git a/src/pspm_sf_dcm.m b/src/pspm_sf_dcm.m
index 1a33f1fea..1fc5b874c 100644
--- a/src/pspm_sf_dcm.m
+++ b/src/pspm_sf_dcm.m
@@ -142,11 +142,17 @@
 end
 miss_epoch = [ymissing_start(:),ymissing_end(:)];
 flag_missing_too_long = 0;
-if any(diff(miss_epoch, 1, 2)/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)/sr > 0)
+  if any(diff(miss_epoch, 1, 2)/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

From f3b5fc566939273f4d07d33dc8fc1dbf844f38dd Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Mon, 29 May 2023 15:35:59 +0100
Subject: [PATCH 02/48] attempting to resolve issue

By changing only internal PsPM code
---
 src/pspm_sf_dcm.m | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/pspm_sf_dcm.m b/src/pspm_sf_dcm.m
index 1fc5b874c..d74ae21ed 100644
--- a/src/pspm_sf_dcm.m
+++ b/src/pspm_sf_dcm.m
@@ -104,12 +104,13 @@
 y = 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));
 priors.muX0 = X0;
-nresp = floor(fresp * numel(y)/sr) + 1;
+nresp = floor(fresp * numel(y_non_nan)/sr) + 1;
 u = [];
 u(1, :) = (1:numel(y))/sr;
 u(2, :) = nresp;

From c0a4a125ba101bba1bfe6aa6eccd1841ea871536 Mon Sep 17 00:00:00 2001
From: "teddy.chao" <dadi.zhao@ucl.ac.uk>
Date: Mon, 5 Jun 2023 14:09:07 +0100
Subject: [PATCH 03/48] update to use non-zero y

---
 src/pspm_dcm_inv.m | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/pspm_dcm_inv.m b/src/pspm_dcm_inv.m
index 7caea6556..08752101f 100644
--- a/src/pspm_dcm_inv.m
+++ b/src/pspm_dcm_inv.m
@@ -551,9 +551,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_nonzero = y(~isnan(y));
+        priors.muX0(1) = mean(y_nonzero(1:3));%) - min(y);
+        priors.muX0(2) = mean(diff(y_nonzero(1:3)));
+        priors.muX0(3) = diff(diff(y_nonzero(1:3)));
         priors.muX0(7) = 0;%min(y);
         for n = [1:3 7]
           priors.SigmaX0(n, n) = 1e-2;

From c1ea66364223e719c069d2d96d5c7d53c5bbec41 Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Mon, 5 Jun 2023 15:41:57 +0100
Subject: [PATCH 04/48] Update pspm_sf_dcm.m

---
 src/pspm_sf_dcm.m | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/pspm_sf_dcm.m b/src/pspm_sf_dcm.m
index d74ae21ed..70c23ad23 100644
--- a/src/pspm_sf_dcm.m
+++ b/src/pspm_sf_dcm.m
@@ -158,7 +158,8 @@
 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

From 0f78326b1de0e87fda89779f9bcd03663a82f54d Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Mon, 5 Jun 2023 16:39:30 +0100
Subject: [PATCH 05/48] Update pspm_dcm_inv.m

---
 src/pspm_dcm_inv.m | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/pspm_dcm_inv.m b/src/pspm_dcm_inv.m
index 65d720f6c..0777d8753 100644
--- a/src/pspm_dcm_inv.m
+++ b/src/pspm_dcm_inv.m
@@ -551,10 +551,10 @@
       priors.SigmaX0 = zeros(7);
       priors.muX0 = zeros(7, 1);
       if trl == 1
-        y_nonzero = y(~isnan(y));
-        priors.muX0(1) = mean(y_nonzero(1:3));%) - min(y);
-        priors.muX0(2) = mean(diff(y_nonzero(1:3)));
-        priors.muX0(3) = diff(diff(y_nonzero(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;

From b677236cddea04cda7e319a3fef46b2af63180c3 Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Mon, 12 Jun 2023 12:55:12 +0100
Subject: [PATCH 06/48] Update pspm_sf_dcm.m

---
 src/pspm_sf_dcm.m | 20 ++++++++++++++++++++
 1 file changed, 20 insertions(+)

diff --git a/src/pspm_sf_dcm.m b/src/pspm_sf_dcm.m
index b2a49460a..b380163ff 100644
--- a/src/pspm_sf_dcm.m
+++ b/src/pspm_sf_dcm.m
@@ -132,6 +132,26 @@
   '\n=========================================================\n'], c(4:6));
 if isfield(options, 'missing')
   ymissing = options.missing;
+  if isfile(ymissing)
+    [~, ~, fExt] = fileparts(ymissing);
+    switch lower(fExt)
+      case '.mat'
+        % A MAT file
+        ymissing = load(ymissing);
+        ymissing = ymissing.epochs;
+      otherwise  % Under all circumstances SWITCH gets an OTHERWISE!
+        warning('ID:invalid_input', 'Missing epochs should be a matlab file if loading from external.');
+    end
+  else
+    if ~isnumeric(ymissing)
+      warning('ID:invalid_input', 'Missing epochs should be numeric if as a variable.');
+    end
+  end
+  if any(size(ymissing)~=size(y))
+    ymiss = ymissing;
+    ymissing = zeros(size(y));
+    ymissing(ymiss(:,1)+1:ymiss(:,2)+1) = 1;
+  end
 else
   ymissing = isnan(y);
 end

From 46d0f0181e66ce16572ed8e25ed10ad699beaddd Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Mon, 12 Jun 2023 14:15:00 +0100
Subject: [PATCH 07/48] Update pspm_split_sessions.m

---
 src/pspm_split_sessions.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pspm_split_sessions.m b/src/pspm_split_sessions.m
index b0c255ffc..4058dc049 100644
--- a/src/pspm_split_sessions.m
+++ b/src/pspm_split_sessions.m
@@ -213,7 +213,7 @@
     pspm_load_data(newdatafile{sn}, newdata);
     % 2.4.5 Split Epochs
     if ~isempty(options.missing)
-      dummydata{1,1}.header = struct('channeltype', 'custom', ...
+      dummydata{1,1}.header = struct('chantype', 'custom', ...
         'sr', missingsr, ...
         'units', 'unknown');
       dummydata{1,1}.data   = dp_epochs;

From 005e97778c7130b46544595c5b124f3af8ccd05b Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Mon, 12 Jun 2023 15:19:29 +0100
Subject: [PATCH 08/48] correct some comment text

---
 src/backroom/pspm_find_data_epochs.m | 18 +++++++++---------
 src/pspm_cfg/pspm_cfg_data_convert.m |  2 +-
 2 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/src/backroom/pspm_find_data_epochs.m b/src/backroom/pspm_find_data_epochs.m
index 2ef36973a..8d7cd0e93 100644
--- a/src/backroom/pspm_find_data_epochs.m
+++ b/src/backroom/pspm_find_data_epochs.m
@@ -3,25 +3,25 @@
 % of a channel. It then writes it to a file depending on the input and
 % returns the location of that file.
 %
-% FORMAT: 
+% FORMAT:
 %   [sts, epochfile] = pspm_find_data_epochs(datafile, chan_id, options)
-% 
+%
 % ARGUMENTS:
 %   datafile:           File which contains the corresponding channel to
 %                       look for epochs.
-%   chan_id:            Number or channeltype of the the channel. If
-%                       multiple channels match the specified channeltype
+%   chan_id:            Number or chantype of the the channel. If
+%                       multiple channels match the specified chantype
 %                       only the first found channel will be used.
-%   options:            
+%   options:
 %       output_file:    Name and path to the output file. Default is the
 %                       same as datafile but prepended with an 'e'.
 %       overwrite:      Overwrite output file if it already exists. Default
 %                       is 0.
-%                       
+%
 % OUTPUT:
 %   sts:                Return status of the function. If equals 1 no
 %                       error or warnings were produced.
-%   epochfile:          File where the epochs have been saved to. 
+%   epochfile:          File where the epochs have been saved to.
 %__________________________________________________________________________
 % PsPM 3.1
 % (C) 2016 Tobias Moser (University of Zurich)
@@ -90,7 +90,7 @@
     if (isempty(stop) && ~isempty(start)) || (start(end) > stop(end))
         stop = [stop; numel(logi)];
     end;
-    
+
     if (isempty(start) && ~isempty(stop)) || (start(1) > stop(1))
         start = [1; start];
     end;
@@ -102,7 +102,7 @@
 file_exist = exist(options.output_file, 'file');
 write_ok = false;
 if file_exist
-    if options.overwrite 
+    if options.overwrite
         write_ok = true;
     elseif ~options.dont_ask_overwrite
         if feature('ShowFigureWindows')
diff --git a/src/pspm_cfg/pspm_cfg_data_convert.m b/src/pspm_cfg/pspm_cfg_data_convert.m
index 428c79fea..62b8d11cd 100644
--- a/src/pspm_cfg/pspm_cfg_data_convert.m
+++ b/src/pspm_cfg/pspm_cfg_data_convert.m
@@ -117,7 +117,7 @@
 visangle2sps.help   = {['Takes paires of channels with gaze data in spherical',...
                        ' coordinates (i.e. visual angle) and computes scanpath speed',...
                        ' (i.e. scalar distance per second). Saves result into a ',...
-                       'new channel with channeltype ''sps'' (scanpath speed)']};
+                       'new channel with chantype ''sps'' (scanpath speed)']};
 
 %% Mode
 mode                = cfg_choice;

From 3ea7e181f0dc6b3439336e7a6edec5121a475ae5 Mon Sep 17 00:00:00 2001
From: "teddy.chao" <dadi.zhao@ucl.ac.uk>
Date: Mon, 19 Jun 2023 11:26:41 +0100
Subject: [PATCH 09/48] update missing allocation

---
 src/pspm_dcm_inv.m | 1 -
 src/pspm_sf.m      | 4 ++++
 2 files changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/pspm_dcm_inv.m b/src/pspm_dcm_inv.m
index a3bdf7908..9c1b47b35 100644
--- a/src/pspm_dcm_inv.m
+++ b/src/pspm_dcm_inv.m
@@ -64,7 +64,6 @@
 %   │                 [numeric, default: 0.1, unit: second]
 %   │                 minimum dispersion (standard deviation) for flexible
 %   │                 responses.
-%   ├──────.missing:  data points to be disregarded by inversion.
 %   │ ▶︎ display
 %   ├──────.dispwin:  [bool, default as 1]
 %   │                 display progress window.
diff --git a/src/pspm_sf.m b/src/pspm_sf.m
index 93fdaadd5..5cc66a71b 100644
--- a/src/pspm_sf.m
+++ b/src/pspm_sf.m
@@ -25,6 +25,7 @@
 %   │                   as 'dcm'
 %   │                   [cell_array] a cell array of methods mentioned above.
 %   ├─────────.filter:  filter settings; modality specific default
+%   ├────────.missing:  index of missing values to ignore
 %   └────────.channel:  channel number; default: first SCR channel
 %   ┌─────────options
 %   ├──────.overwrite:  [logical] (0 or 1)
@@ -170,6 +171,9 @@
 end
 % 2.8 Set options
 try model.channel; catch, model.channel = 'scr'; end
+if isfield(model, 'missing')
+  options.missing = model.missing;
+end
 options = pspm_options(options, 'sf');
 if options.invalid
   return

From 9da282caec883c17a6b3efdcdb12f70ab389cb9c Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Thu, 22 Jun 2023 15:02:32 +0100
Subject: [PATCH 10/48] Update pspm_sf_dcm.m

---
 src/pspm_sf_dcm.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pspm_sf_dcm.m b/src/pspm_sf_dcm.m
index b380163ff..1971ccd44 100644
--- a/src/pspm_sf_dcm.m
+++ b/src/pspm_sf_dcm.m
@@ -111,7 +111,7 @@
 X0(2, 1)   = mean(diff(x0));
 X0(3, 1)   = diff(diff(x0));
 priors.muX0 = X0;
-nresp = floor(fresp * numel(y_non_nan)/sr) + 1;
+nresp = floor(fresp * numel(y)/sr) + 1;
 u = [];
 u(1, :) = (1:numel(y))/sr;
 u(2, :) = nresp;

From 68710990291a45818ed6f2fa8e3bd1affc544ea7 Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Thu, 22 Jun 2023 15:44:11 +0100
Subject: [PATCH 11/48] Update pspm_dcm_inv.m

---
 src/pspm_dcm_inv.m | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/pspm_dcm_inv.m b/src/pspm_dcm_inv.m
index 9c1b47b35..a03b5a1d0 100644
--- a/src/pspm_dcm_inv.m
+++ b/src/pspm_dcm_inv.m
@@ -34,6 +34,7 @@
 %   │                 transformed back into raw data units.
 %   ├───.flexevents:  flexible events to adjust amplitude priors
 %   ├────.fixevents:  fixed events to adjust amplitude priors
+%   ├─.missing_data:  file storing missing epochs
 %   └──.constrained:  [optional]
 %                     constrained model for flexible responses which have
 %                     have fixed dispersion (0.3 s SD) but flexible latency
@@ -125,9 +126,9 @@
 try model.meanSCR; catch, model.meanSCR = 0; end
 % These parameters were set with default fallback values but will be
 % determined later by processing (same to pspm_dcm)
-%try model.fixevents; catch, warning('model.fixevents not defined.'); end
-%try model.flexevents; catch, warning('model.flexevents not defined.'); end
-%try model.missing_data; catch, warning('model.missing_data not defined.'); end
+try model.fixevents; catch, warning('model.fixevents not defined.'); end
+try model.flexevents; catch, warning('model.flexevents not defined.'); end
+try model.missing_data; catch, warning('model.missing_data not defined.'); end
 % These parameters do not need to have a default value and will be
 % determined later (same to pspm_dcm)
 

From 6dfe537d4fdfc26d6674d713c18b04022f504996 Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Thu, 22 Jun 2023 21:44:26 +0100
Subject: [PATCH 12/48] Update pspm_sf_dcm.m

---
 src/pspm_sf_dcm.m | 41 ++++++++++++++++++++++++++---------------
 1 file changed, 26 insertions(+), 15 deletions(-)

diff --git a/src/pspm_sf_dcm.m b/src/pspm_sf_dcm.m
index 494a854f1..fb4a8d3fd 100644
--- a/src/pspm_sf_dcm.m
+++ b/src/pspm_sf_dcm.m
@@ -1,4 +1,4 @@
-function varargout = pspm_sf_dcm(scr, sr, options)
+function varargout = pspm_sf_dcm(model, options)
 % ● Description
 %   pspm_sf_dcm does dynamic causal modelling for SF of the skin conductance
 %   uses f_SF and g_Id
@@ -18,19 +18,24 @@
 %       .yhat:  fitted time series
 %      .model:  information about the DCM inversion
 % ● Arguments
-%         scr:  skin conductance epoch (maximum size depends on computing
-%               power, a sensible size is 60 s at 10 Hz)
-%          sr:  sampling rate in Hz
-%     options:  options structure
-%  .threshold:  threshold for SN detection (default 0.1 mcS)
-%      .theta:  a (1 x 5) vector of theta values for f_SF
-%               (default: read from pspm_sf_theta)
-%      .fresp:  frequency of responses to model (default 0.5 Hz)
-%    .dispwin:  display progress window (default 1)
-%   .dispsmallwin:
-%               display intermediate windows (default 0);
-%    .missing:  index of missing values to ignore
-%   .missingthresh:
+%   ┌──────model
+%   │ ▶︎ Mandatory
+%   ├───────.scr:  skin conductance epoch (maximum size depends on computing
+%   │              power, a sensible size is 60 s at 10 Hz)
+%   ├────────.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
+%   │              (default: read from pspm_sf_theta)
+%   ├─────.fresp:  frequency of responses to model (default 0.5 Hz)
+%   ├───.dispwin:  display progress window (default 1)
+%   ├.dispsmallwin:
+%   │              display intermediate windows (default 0);
+%   ├───.missing:  index of missing values to ignore
+%   └─.missingthresh:
 %               threshold value for controlling missing epochs (default 2s).
 % ● References
 %   Bach DR, Daunizeau J, Kuelzow N, Friston KJ, & Dolan RJ (2011). Dynamic
@@ -47,8 +52,14 @@
 end
 sts = -1;
 tstart = tic;
+
 %% 2 Check input arguments
-if nargin < 2 || ~isnumeric(sr) || numel(sr) > 1
+% 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(sr) || numel(sr) > 1
   errmsg = sprintf('No valid sample rate given.');
 elseif (sr < 1) || (sr > 1e5)
   errmsg = sprintf('Sample rate out of range.');

From a28018bab9b97c4af56d9ef0dcba3090839c86fb Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Thu, 22 Jun 2023 22:15:00 +0100
Subject: [PATCH 13/48] update sf and sf_dcm

---
 src/pspm_sf.m     |  3 ++-
 src/pspm_sf_dcm.m | 34 +++++++++++++++++-----------------
 2 files changed, 19 insertions(+), 18 deletions(-)

diff --git a/src/pspm_sf.m b/src/pspm_sf.m
index 5cc66a71b..02a8e0235 100644
--- a/src/pspm_sf.m
+++ b/src/pspm_sf.m
@@ -263,7 +263,8 @@
       escr = Y{datatype(k)}(win(1):win(end));
       sf.model{k}(iEpoch).data = escr;
       % 3.7 do the analysis and collect results
-      invrs = fhandle{k}(escr, sr(datatype(k)), options);
+      model_analysis = struct('scr', escr, 'sr', sr(datatype(k)));
+      invrs = fhandle{k}(model_analysis, options);
       if any(strcmpi(method{k}, {'dcm', 'mp'}))
         sf.model{k}(iEpoch).inv     = invrs;
         sf.stats(iEpoch, k)         = invrs.f;
diff --git a/src/pspm_sf_dcm.m b/src/pspm_sf_dcm.m
index fb4a8d3fd..743c3cec0 100644
--- a/src/pspm_sf_dcm.m
+++ b/src/pspm_sf_dcm.m
@@ -4,7 +4,7 @@
 %   uses f_SF and g_Id
 %   the input data is assumed to be in mcS, and sampling rate in Hz
 % ● Format
-%   function out = pspm_sf_dcm(scr, sr, options)
+%   function out = pspm_sf_dcm(model, options)
 % ● Output
 %         out:  output
 %          .n:  number of responses above threshold
@@ -59,18 +59,18 @@
 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(sr) || numel(sr) > 1
+if ~isnumeric(model.sr) || numel(model.sr) > 1
   errmsg = sprintf('No valid sample rate given.');
-elseif (sr < 1) || (sr > 1e5)
+elseif (model.sr < 1) || (model.sr > 1e5)
   errmsg = sprintf('Sample rate out of range.');
-elseif exist('osr', 'var') && osr ~= sr
+elseif exist('osr', 'var') && osr ~= model.sr
   errmsg = sprintf('Sample rate of theta file is different from sample rate of data.');
-elseif nargin < 1 || ~isnumeric(scr)
+elseif nargin < 1 || ~isnumeric(model.scr)
   errmsg = 'No data.';
-elseif ~any(size(scr) == 1)
+elseif ~any(size(model.scr) == 1)
   errmsg = 'Input SCR is not a vector';
 else
-  scr = scr(:);
+  model.scr = model.scr(:);
 end
 if exist('errmsg', 'var') == 1
   warning(errmsg);
@@ -104,16 +104,16 @@
 priors.a_alpha = Inf;
 priors.b_alpha = 0;
 % 4.2 initialise priors in correct dimensions
-priors.iQy = cell(numel(scr), 1);
-priors.iQx = cell(numel(scr), 1);
-for k = 1:numel(scr)  % default priors on noise covariance
+priors.iQy = cell(numel(model.scr), 1);
+priors.iQx = cell(numel(model.scr), 1);
+for k = 1:numel(model.scr)  % default priors on noise covariance
   priors.iQy{k} = 1;
   priors.iQx{k} = eye(dim.n);
 end
 options.inG.ind = 1;
-options.inF.dt = 1/sr;
+options.inF.dt = 1/model.sr;
 % 4.3 prepare data
-y = scr;
+y = model.scr;
 y = y - min(y);
 % 4.4 determine initial conditions
 x0 = y(1:3);
@@ -121,9 +121,9 @@
 X0(2, 1)   = mean(diff(x0));
 X0(3, 1)   = diff(diff(x0));
 priors.muX0 = X0;
-nresp = floor(fresp * numel(y)/sr) + 1;
+nresp = floor(fresp * numel(y)/model.sr) + 1;
 u = [];
-u(1, :) = (1:numel(y))/sr;
+u(1, :) = (1:numel(y))/model.sr;
 u(2, :) = nresp;
 priors.muTheta = theta(1:3)';
 priors.muTheta(4:2:(2 * nresp + 3)) = 1/fresp * (0:(nresp-1));
@@ -154,7 +154,7 @@
 end
 miss_epoch = [ymissing_start(:),ymissing_end(:)];
 flag_missing_too_long = 0;
-if any(diff(miss_epoch, 1, 2)/sr > options.missingthresh)
+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('ID:missing_data', warning_message);
@@ -169,7 +169,7 @@
   end
   t = posterior.muTheta(4:2:end);
   a = exp(posterior.muTheta(5:2:end) - theta(5));   % rescale
-  ex = find(t < -2 | t > (numel(scr)/sr - 1)); % find SA responses the SCR peak of which is outside episode
+  ex = find(t < -2 | t > (numel(model.scr)/model.sr - 1)); % find SA responses the SCR peak of which is outside episode
   t(ex) = [];
   a(ex) = [];
 end
@@ -178,7 +178,7 @@
   out.t               = t - theta(4);   % subtract conduction delay
   out.a               = a;
   out.n               = numel(find(a > threshold));
-  out.f               = out.n/(numel(scr)/sr);
+  out.f               = out.n/(numel(model.scr)/model.sr);
   out.ma              = mean(a(a > threshold));
   out.theta           = theta;
   out.if              = fresp;

From 8b34c7e376e9499287225c55fadad73c58988598 Mon Sep 17 00:00:00 2001
From: "teddy.chao" <dadi.zhao@ucl.ac.uk>
Date: Sat, 24 Jun 2023 19:30:03 +0100
Subject: [PATCH 14/48] update missing and missing_data description

---
 src/pspm_dcm.m     | 2 +-
 src/pspm_dcm_inv.m | 4 +++-
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/pspm_dcm.m b/src/pspm_dcm.m
index d668beb8e..732368585 100644
--- a/src/pspm_dcm.m
+++ b/src/pspm_dcm.m
@@ -30,7 +30,7 @@
 %   │             trial. If this is not the case, include `dummy` events with
 %   │             negative onsets.
 %   │ ▶︎ Optional
-%   ├───.missing: Allows to specify missing (e. g. artefact) epochs in the
+%   ├───.missing: Allows to specify missing (e.g. artefact) epochs in the
 %   │             data file. See pspm_get_timing for epoch definition; specify
 %   │             a cell array for multiple input files. This must always be
 %   │             specified in SECONDS.
diff --git a/src/pspm_dcm_inv.m b/src/pspm_dcm_inv.m
index a03b5a1d0..30b130204 100644
--- a/src/pspm_dcm_inv.m
+++ b/src/pspm_dcm_inv.m
@@ -34,7 +34,9 @@
 %   │                 transformed back into raw data units.
 %   ├───.flexevents:  flexible events to adjust amplitude priors
 %   ├────.fixevents:  fixed events to adjust amplitude priors
-%   ├─.missing_data:  file storing missing epochs
+%   ├─.missing_data:  missing epoch data, originally loaded as model.missing
+%   │                 from pspm_dcm, but calculated into .missing_data in pspm_dcm
+%   │                 and then transferred to pspm_dcm_inv.
 %   └──.constrained:  [optional]
 %                     constrained model for flexible responses which have
 %                     have fixed dispersion (0.3 s SD) but flexible latency

From 8a8a77a7d895ac9f5d6472b277c021a696145915 Mon Sep 17 00:00:00 2001
From: "teddy.chao" <dadi.zhao@ucl.ac.uk>
Date: Sat, 24 Jun 2023 21:13:34 +0100
Subject: [PATCH 15/48] update missing checking in sf

---
 src/pspm_dcm_inv.m |  4 ++--
 src/pspm_sf.m      | 45 ++++++++++++++++++++++++++++++++++++++-------
 2 files changed, 40 insertions(+), 9 deletions(-)

diff --git a/src/pspm_dcm_inv.m b/src/pspm_dcm_inv.m
index 30b130204..7066e8ad7 100644
--- a/src/pspm_dcm_inv.m
+++ b/src/pspm_dcm_inv.m
@@ -35,8 +35,8 @@
 %   ├───.flexevents:  flexible events to adjust amplitude priors
 %   ├────.fixevents:  fixed events to adjust amplitude priors
 %   ├─.missing_data:  missing epoch data, originally loaded as model.missing
-%   │                 from pspm_dcm, but calculated into .missing_data in pspm_dcm
-%   │                 and then transferred to pspm_dcm_inv.
+%   │                 from pspm_dcm, but calculated into .missing_data (created
+%   │                 in pspm_dcm and then transferred to pspm_dcm_inv.
 %   └──.constrained:  [optional]
 %                     constrained model for flexible responses which have
 %                     have fixed dispersion (0.3 s SD) but flexible latency
diff --git a/src/pspm_sf.m b/src/pspm_sf.m
index 02a8e0235..384b6d4cf 100644
--- a/src/pspm_sf.m
+++ b/src/pspm_sf.m
@@ -25,7 +25,12 @@
 %   │                   as 'dcm'
 %   │                   [cell_array] a cell array of methods mentioned above.
 %   ├─────────.filter:  filter settings; modality specific default
-%   ├────────.missing:  index of missing values to ignore
+%   ├────────.missing:  [string/cell_array]
+%   │                   Allows to specify missing (e.g. artefact) epochs in the
+%   │                   data file. See pspm_get_timing for epoch definition; specify
+%   │                   a cell array for multiple input files. This must always be
+%   │                   specified in SECONDS.
+%   │                   Default: no missing values
 %   └────────.channel:  channel number; default: first SCR channel
 %   ┌─────────options
 %   ├──────.overwrite:  [logical] (0 or 1)
@@ -171,13 +176,20 @@
 end
 % 2.8 Set options
 try model.channel; catch, model.channel = 'scr'; end
-if isfield(model, 'missing')
-  options.missing = model.missing;
-end
 options = pspm_options(options, 'sf');
 if options.invalid
   return
 end
+% 2.9 Set missing epochs
+if ~isfield(model, 'missing')
+  model.missing = cell(numel(model.datafile), 1);
+elseif ischar(model.missing) || isnumeric(model.missing)
+  model.missing = {model.missing};
+elseif ~iscell(model.missing)
+  warning('ID:invalid_input',...
+    'Missing values must be a filename, matrix, or cell array of these.');
+  return
+end
 %% 3 Get data
 for iFile = 1:numel(model.datafile)
   % 3.1 User output
@@ -203,7 +215,26 @@
       'Please check your results.\n'], ...
       data{1}.header.units);
   end
-  % 3.5 Get marker data
+  % 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,:);
+      % 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, :) = [];
+        end
+      end
+    end
+  else
+    missing{iFile} = [];
+  end
+  % 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);
@@ -253,7 +284,7 @@
         win(1) = max(win(1), 1);
         win(2) = min(win(2), numel(Y{datatype(k)}));
       end
-      % 3.6 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;
@@ -262,7 +293,7 @@
       %
       escr = Y{datatype(k)}(win(1):win(end));
       sf.model{k}(iEpoch).data = escr;
-      % 3.7 do the analysis and collect results
+      % 3.6.2 do the analysis and collect results
       model_analysis = struct('scr', escr, 'sr', sr(datatype(k)));
       invrs = fhandle{k}(model_analysis, options);
       if any(strcmpi(method{k}, {'dcm', 'mp'}))

From 4200b2a716e2ec548bbc3d6225a417c0a07d5915 Mon Sep 17 00:00:00 2001
From: "teddy.chao" <dadi.zhao@ucl.ac.uk>
Date: Sun, 25 Jun 2023 18:16:23 +0100
Subject: [PATCH 16/48] Update pspm_sf_dcm.m

---
 src/pspm_sf_dcm.m | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/pspm_sf_dcm.m b/src/pspm_sf_dcm.m
index aed4ba771..36ca689c7 100644
--- a/src/pspm_sf_dcm.m
+++ b/src/pspm_sf_dcm.m
@@ -175,8 +175,8 @@
 end
 miss_epoch = [ymissing_start(:),ymissing_end(:)];
 flag_missing_too_long = 0;
-if any(diff(miss_epoch, 1, 2)/sr > 0)
-  if any(diff(miss_epoch, 1, 2)/sr > options.missingthresh)
+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;

From 9b5cb8573559c86b0c6beb12c000ce8b0bf81332 Mon Sep 17 00:00:00 2001
From: "teddy.chao" <dadi.zhao@ucl.ac.uk>
Date: Sun, 25 Jun 2023 18:59:09 +0100
Subject: [PATCH 17/48] update missing field reading

---
 src/pspm_dcm_inv.m | 6 +++---
 src/pspm_sf_dcm.m  | 5 ++---
 2 files changed, 5 insertions(+), 6 deletions(-)

diff --git a/src/pspm_dcm_inv.m b/src/pspm_dcm_inv.m
index 7066e8ad7..51d9df16a 100644
--- a/src/pspm_dcm_inv.m
+++ b/src/pspm_dcm_inv.m
@@ -128,9 +128,9 @@
 try model.meanSCR; catch, model.meanSCR = 0; end
 % These parameters were set with default fallback values but will be
 % determined later by processing (same to pspm_dcm)
-try model.fixevents; catch, warning('model.fixevents not defined.'); end
-try model.flexevents; catch, warning('model.flexevents not defined.'); end
-try model.missing_data; catch, warning('model.missing_data not defined.'); end
+% try model.fixevents; catch, warning('model.fixevents not defined.'); end
+% try model.flexevents; catch, warning('model.flexevents not defined.'); end
+% try model.missing_data; catch, warning('model.missing_data not defined.'); end
 % These parameters do not need to have a default value and will be
 % determined later (same to pspm_dcm)
 
diff --git a/src/pspm_sf_dcm.m b/src/pspm_sf_dcm.m
index 743c3cec0..688b59ba8 100644
--- a/src/pspm_sf_dcm.m
+++ b/src/pspm_sf_dcm.m
@@ -34,7 +34,6 @@
 %   ├───.dispwin:  display progress window (default 1)
 %   ├.dispsmallwin:
 %   │              display intermediate windows (default 0);
-%   ├───.missing:  index of missing values to ignore
 %   └─.missingthresh:
 %               threshold value for controlling missing epochs (default 2s).
 % ● References
@@ -140,8 +139,8 @@
 c = clock;
 fprintf(['\n\nEstimating model parameters for f_SF ... \t%02.0f:%02.0f:%02.0f', ...
   '\n=========================================================\n'], c(4:6));
-if isfield(options, 'missing')
-  ymissing = options.missing;
+if isfield(model, 'missing_data')
+  ymissing = model.missing_data;
 else
   ymissing = isnan(y);
 end

From 4fde373065d775c8a0332a6276d006acf2f84cd7 Mon Sep 17 00:00:00 2001
From: "teddy.chao" <dadi.zhao@ucl.ac.uk>
Date: Sun, 25 Jun 2023 23:53:21 +0100
Subject: [PATCH 18/48] update SF

---
 src/pspm_sf.m     | 126 +++++++++++++++++++++++++++++++++++++---------
 src/pspm_sf_dcm.m |  31 ++++++------
 2 files changed, 115 insertions(+), 42 deletions(-)

diff --git a/src/pspm_sf.m b/src/pspm_sf.m
index 384b6d4cf..12102065b 100644
--- a/src/pspm_sf.m
+++ b/src/pspm_sf.m
@@ -191,15 +191,15 @@
   return
 end
 %% 3 Get data
-for iFile = 1:numel(model.datafile)
+for iSn = 1:numel(model.datafile)
   % 3.1 User output
-  fprintf('SF analysis: %s ...', model.datafile{iFile});
+  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);
+  [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,31 +208,107 @@
     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.5.3 find missing epochs according to subsession threshold ---
+  n_data = size(data{iSn}{1}.data,1);
+  if isempty(missing{iSn})
+    nan_epochs = isnan(data{iSn}{1}.data);
+    d_nan_ep = transpose(diff(nan_epochs));
+    nan_ep_start = find(d_nan_ep == 1);
+    nan_ep_stop = find(d_nan_ep == -1);
+    if numel(nan_ep_start) > 0 || numel(nan_ep_stop) > 0
+      % check for blunt ends and fix
+      if isempty(nan_ep_start)
+        nan_ep_start = 1;
+      elseif isempty(nan_ep_stop)
+        nan_ep_stop = numel(d_nan_ep);
+      end
+      if nan_ep_start(1) > nan_ep_stop(1)
+        nan_ep_start = [1, nan_ep_start];
+      end
+      if nan_ep_start(end) > nan_ep_stop(end)
+        nan_ep_stop(end + 1) = numel(d_nan_ep);
+      end
+    end
+    % put missing epochs together
+    miss_epochs = [nan_ep_start(:), nan_ep_stop(:)];
+    % classify if epoch should be considered
+    % true for duration > substhresh and for missing epochs
+    ignore_epochs = diff(miss_epochs, 1, 2)/data{iSn}{1}.header.sr > ...
+      model.substhresh;
+  else
+    % use missing epochs as specified by file
+    miss_epochs = pspm_time2index(missing{iSn},data{iSn}{1}.header.sr);
+    ignore_epochs = diff(missing{iSn}, 1, 2) > model.substhresh;
+    % and set data to NaN to enable later detection of `short` missing
+    % epochs
+    for k = 1:size(miss_epochs, 1)
+      flanks = round(miss_epochs(k,:));
+      data{iSn}{1}.data(flanks(1):flanks(2)) = NaN;
+    end
+  end
+  if any(ignore_epochs)
+    i_e = find(ignore_epochs);
+    % invert missings to sessions without nans
+    se_start = [1; miss_epochs(i_e(1:end), 2) + 1];
+    se_stop = [miss_epochs(i_e(1:end), 1)-1; n_data];
+    % throw away first session if stop is
+    % earlier than start (can happen because stop - 1)
+    % is used
+    if se_stop(1) <= se_start(1)
+      se_start = se_start(2:end);
+      se_stop = se_stop(2:end);
+    end
+    % throw away last session if start (+1) overlaps
+    % n_data
+    if se_start(end) >= n_data
+      se_start = se_start(1:end-1);
+      se_stop = se_stop(1:end-1);
+    end
+    % subsessions header --
+    % =====================
+    % 1 session_id
+    % 2 start_time (s)
+    % 3 stop_time (s)
+    % 4 missing (1) or data segment (0)
+    n_sbs = numel(se_start);
+    % enabled subsessions
+    subsessions(end+(1:n_sbs), 1:4) = [ones(n_sbs,1)*iSn, ...
+      [se_start, se_stop]/data{iSn}{1}.header.sr, ...
+      zeros(n_sbs,1)];
+    % missing epochs
+    n_miss = sum(ignore_epochs);
+    subsessions(end+(1:n_miss), 1:4) = [ones(n_miss,1)*iSn, ...
+      miss_epochs(i_e,:)/data{iSn}{1}.header.sr, ...
+      ones(n_miss,1)];
+  else
+    subsessions(end+1,1:4) = [iSn, ...
+      [1, numel(data{iSn}{1}.data)]/data{iSn}{1}.header.sr, 0];
   end
   % 3.6 Get marker data
   if any(strcmp(model.timeunits, {'marker', 'markers'}))
@@ -262,23 +338,23 @@
     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);
@@ -286,7 +362,7 @@
       end
       % 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));
@@ -308,16 +384,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;
diff --git a/src/pspm_sf_dcm.m b/src/pspm_sf_dcm.m
index 36ca689c7..104cf92a7 100644
--- a/src/pspm_sf_dcm.m
+++ b/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
@@ -37,6 +24,18 @@
 %   ├───.missing:  index of missing values to ignore
 %   └─.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.
@@ -52,12 +51,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.');
@@ -126,14 +123,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;

From 0702106698e0cca1d88797042a34542a22d52a43 Mon Sep 17 00:00:00 2001
From: "teddy.chao" <dadi.zhao@ucl.ac.uk>
Date: Mon, 26 Jun 2023 01:20:26 +0100
Subject: [PATCH 19/48] Update pspm_sf.m

---
 src/pspm_sf.m | 52 +++++++++++++++++++++++++++++----------------------
 1 file changed, 30 insertions(+), 22 deletions(-)

diff --git a/src/pspm_sf.m b/src/pspm_sf.m
index 12102065b..96712db11 100644
--- a/src/pspm_sf.m
+++ b/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,7 +109,7 @@
 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;
@@ -117,7 +117,14 @@
   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 substhresh --
+if ~isfield(model, 'substhresh')
+  model.substhresh = 2;
+elseif ~isnumeric(model.substhresh)
+  warning('ID:invalid_input', 'Subsession threshold must be numeric.');
+  return;
+end
+% 2.6 check methods --
 if ~isfield(model, 'method')
   model.method = {'dcm'};
 elseif ischar(model.method)
@@ -156,7 +163,7 @@
     end
   end
 end
-% 2.6 Check timing
+% 2.7 Check timing --
 if strcmpi(model.timeunits, 'whole')
   epochs = repmat({[1 1]}, numel(model.datafile), 1);
 else
@@ -168,19 +175,19 @@
     end
   end
 end
-% 2.7 Check filter
+% 2.8 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.9 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.10 Set missing epochs --
 if ~isfield(model, 'missing')
   model.missing = cell(numel(model.datafile), 1);
 elseif ischar(model.missing) || isnumeric(model.missing)
@@ -191,14 +198,15 @@
   return
 end
 %% 3 Get data
+subsessions = [];
 for iSn = 1:numel(model.datafile)
-  % 3.1 User output
+  % 3.1 User output --
   fprintf('SF analysis: %s ...', model.datafile{iSn});
-  % 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{iSn}, model.channel);
   if sts_load_data < 0, return; end
   Y{1} = data{1}.data; sr(1) = data{1}.header.sr;
@@ -208,14 +216,14 @@
     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{iSn})
     [~, missing{iSn}] = pspm_get_timing('epochs', model.missing{iSn}, 'seconds');
@@ -235,9 +243,9 @@
     missing{iSn} = [];
   end
   % 3.5.3 find missing epochs according to subsession threshold ---
-  n_data = size(data{iSn}{1}.data,1);
+  n_data = size(data{iSn}.data,1);
   if isempty(missing{iSn})
-    nan_epochs = isnan(data{iSn}{1}.data);
+    nan_epochs = isnan(data{iSn}.data);
     d_nan_ep = transpose(diff(nan_epochs));
     nan_ep_start = find(d_nan_ep == 1);
     nan_ep_stop = find(d_nan_ep == -1);
@@ -259,17 +267,17 @@
     miss_epochs = [nan_ep_start(:), nan_ep_stop(:)];
     % classify if epoch should be considered
     % true for duration > substhresh and for missing epochs
-    ignore_epochs = diff(miss_epochs, 1, 2)/data{iSn}{1}.header.sr > ...
+    ignore_epochs = diff(miss_epochs, 1, 2)/data{iSn}.header.sr > ...
       model.substhresh;
   else
     % use missing epochs as specified by file
-    miss_epochs = pspm_time2index(missing{iSn},data{iSn}{1}.header.sr);
+    miss_epochs = pspm_time2index(missing{iSn},data{iSn}.header.sr);
     ignore_epochs = diff(missing{iSn}, 1, 2) > model.substhresh;
     % and set data to NaN to enable later detection of `short` missing
     % epochs
     for k = 1:size(miss_epochs, 1)
       flanks = round(miss_epochs(k,:));
-      data{iSn}{1}.data(flanks(1):flanks(2)) = NaN;
+      data{iSn}.data(flanks(1):flanks(2)) = NaN;
     end
   end
   if any(ignore_epochs)
@@ -290,7 +298,7 @@
       se_start = se_start(1:end-1);
       se_stop = se_stop(1:end-1);
     end
-    % subsessions header --
+    % subsessions header 
     % =====================
     % 1 session_id
     % 2 start_time (s)
@@ -299,16 +307,16 @@
     n_sbs = numel(se_start);
     % enabled subsessions
     subsessions(end+(1:n_sbs), 1:4) = [ones(n_sbs,1)*iSn, ...
-      [se_start, se_stop]/data{iSn}{1}.header.sr, ...
+      [se_start, se_stop]/data{iSn}.header.sr, ...
       zeros(n_sbs,1)];
     % missing epochs
     n_miss = sum(ignore_epochs);
     subsessions(end+(1:n_miss), 1:4) = [ones(n_miss,1)*iSn, ...
-      miss_epochs(i_e,:)/data{iSn}{1}.header.sr, ...
+      miss_epochs(i_e,:)/data{iSn}.header.sr, ...
       ones(n_miss,1)];
   else
     subsessions(end+1,1:4) = [iSn, ...
-      [1, numel(data{iSn}{1}.data)]/data{iSn}{1}.header.sr, 0];
+      [1, numel(data{iSn}.data)]/data{iSn}.header.sr, 0];
   end
   % 3.6 Get marker data
   if any(strcmp(model.timeunits, {'marker', 'markers'}))

From 3b7f8feb6c592711de29d3651d3417ccd998d9bf Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Mon, 26 Jun 2023 02:05:48 +0100
Subject: [PATCH 20/48] load missing to sf_dcm

---
 src/pspm_sf.m     | 82 +++--------------------------------------------
 src/pspm_sf_dcm.m | 20 ------------
 2 files changed, 4 insertions(+), 98 deletions(-)

diff --git a/src/pspm_sf.m b/src/pspm_sf.m
index 96712db11..e1ea16492 100644
--- a/src/pspm_sf.m
+++ b/src/pspm_sf.m
@@ -198,7 +198,7 @@
   return
 end
 %% 3 Get data
-subsessions = [];
+missing = cell(size(model.missing));
 for iSn = 1:numel(model.datafile)
   % 3.1 User output --
   fprintf('SF analysis: %s ...', model.datafile{iSn});
@@ -242,82 +242,6 @@
   else
     missing{iSn} = [];
   end
-  % 3.5.3 find missing epochs according to subsession threshold ---
-  n_data = size(data{iSn}.data,1);
-  if isempty(missing{iSn})
-    nan_epochs = isnan(data{iSn}.data);
-    d_nan_ep = transpose(diff(nan_epochs));
-    nan_ep_start = find(d_nan_ep == 1);
-    nan_ep_stop = find(d_nan_ep == -1);
-    if numel(nan_ep_start) > 0 || numel(nan_ep_stop) > 0
-      % check for blunt ends and fix
-      if isempty(nan_ep_start)
-        nan_ep_start = 1;
-      elseif isempty(nan_ep_stop)
-        nan_ep_stop = numel(d_nan_ep);
-      end
-      if nan_ep_start(1) > nan_ep_stop(1)
-        nan_ep_start = [1, nan_ep_start];
-      end
-      if nan_ep_start(end) > nan_ep_stop(end)
-        nan_ep_stop(end + 1) = numel(d_nan_ep);
-      end
-    end
-    % put missing epochs together
-    miss_epochs = [nan_ep_start(:), nan_ep_stop(:)];
-    % classify if epoch should be considered
-    % true for duration > substhresh and for missing epochs
-    ignore_epochs = diff(miss_epochs, 1, 2)/data{iSn}.header.sr > ...
-      model.substhresh;
-  else
-    % use missing epochs as specified by file
-    miss_epochs = pspm_time2index(missing{iSn},data{iSn}.header.sr);
-    ignore_epochs = diff(missing{iSn}, 1, 2) > model.substhresh;
-    % and set data to NaN to enable later detection of `short` missing
-    % epochs
-    for k = 1:size(miss_epochs, 1)
-      flanks = round(miss_epochs(k,:));
-      data{iSn}.data(flanks(1):flanks(2)) = NaN;
-    end
-  end
-  if any(ignore_epochs)
-    i_e = find(ignore_epochs);
-    % invert missings to sessions without nans
-    se_start = [1; miss_epochs(i_e(1:end), 2) + 1];
-    se_stop = [miss_epochs(i_e(1:end), 1)-1; n_data];
-    % throw away first session if stop is
-    % earlier than start (can happen because stop - 1)
-    % is used
-    if se_stop(1) <= se_start(1)
-      se_start = se_start(2:end);
-      se_stop = se_stop(2:end);
-    end
-    % throw away last session if start (+1) overlaps
-    % n_data
-    if se_start(end) >= n_data
-      se_start = se_start(1:end-1);
-      se_stop = se_stop(1:end-1);
-    end
-    % subsessions header 
-    % =====================
-    % 1 session_id
-    % 2 start_time (s)
-    % 3 stop_time (s)
-    % 4 missing (1) or data segment (0)
-    n_sbs = numel(se_start);
-    % enabled subsessions
-    subsessions(end+(1:n_sbs), 1:4) = [ones(n_sbs,1)*iSn, ...
-      [se_start, se_stop]/data{iSn}.header.sr, ...
-      zeros(n_sbs,1)];
-    % missing epochs
-    n_miss = sum(ignore_epochs);
-    subsessions(end+(1:n_miss), 1:4) = [ones(n_miss,1)*iSn, ...
-      miss_epochs(i_e,:)/data{iSn}.header.sr, ...
-      ones(n_miss,1)];
-  else
-    subsessions(end+1,1:4) = [iSn, ...
-      [1, numel(data{iSn}.data)]/data{iSn}.header.sr, 0];
-  end
   % 3.6 Get marker data
   if any(strcmp(model.timeunits, {'marker', 'markers'}))
     if options.marker_chan_num
@@ -377,8 +301,10 @@
       %
       escr = Y{datatype(k)}(win(1):win(end));
       sf.model{k}(iEpoch).data = escr;
+      model.missing_data = zeros(size(escr));
+      model.missing_data((missing{iSn}(:,1)+1):(missing{iSn}(:,2)+1)) = 1;
       % 3.6.2 do the analysis and collect results
-      model_analysis = struct('scr', escr, 'sr', sr(datatype(k)));
+      model_analysis = struct('scr', escr, 'sr', sr(datatype(k)), 'missing_data', model.missing_data);
       invrs = fhandle{k}(model_analysis, options);
       if any(strcmpi(method{k}, {'dcm', 'mp'}))
         sf.model{k}(iEpoch).inv     = invrs;
diff --git a/src/pspm_sf_dcm.m b/src/pspm_sf_dcm.m
index 30b7d9c85..482673288 100644
--- a/src/pspm_sf_dcm.m
+++ b/src/pspm_sf_dcm.m
@@ -139,26 +139,6 @@
   '\n=========================================================\n'], c(4:6));
 if isfield(model, 'missing_data')
   ymissing = model.missing_data;
-  if isfile(ymissing)
-    [~, ~, fExt] = fileparts(ymissing);
-    switch lower(fExt)
-      case '.mat'
-        % A MAT file
-        ymissing = load(ymissing);
-        ymissing = ymissing.epochs;
-      otherwise  % Under all circumstances SWITCH gets an OTHERWISE!
-        warning('ID:invalid_input', 'Missing epochs should be a matlab file if loading from external.');
-    end
-  else
-    if ~isnumeric(ymissing)
-      warning('ID:invalid_input', 'Missing epochs should be numeric if as a variable.');
-    end
-  end
-  if any(size(ymissing)~=size(y))
-    ymiss = ymissing;
-    ymissing = zeros(size(y));
-    ymissing(ymiss(:,1)+1:ymiss(:,2)+1) = 1;
-  end
 else
   ymissing = isnan(y);
 end

From 4260862684336c800583c14c338cd8659d6ae748 Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Mon, 26 Jun 2023 02:39:25 +0100
Subject: [PATCH 21/48] Update pspm_sf.m

---
 src/pspm_sf.m | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/src/pspm_sf.m b/src/pspm_sf.m
index e1ea16492..0cad00a31 100644
--- a/src/pspm_sf.m
+++ b/src/pspm_sf.m
@@ -301,10 +301,16 @@
       %
       escr = Y{datatype(k)}(win(1):win(end));
       sf.model{k}(iEpoch).data = escr;
-      model.missing_data = zeros(size(escr));
-      model.missing_data((missing{iSn}(:,1)+1):(missing{iSn}(:,2)+1)) = 1;
+      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
-      model_analysis = struct('scr', escr, 'sr', sr(datatype(k)), 'missing_data', model.missing_data);
+      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;

From 7537b6d97895fe6ff36d7423562d40ab5a5926cd Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Mon, 26 Jun 2023 12:35:58 +0100
Subject: [PATCH 22/48] Update pspm_sf.m

---
 src/pspm_sf.m | 27 ++++++++++-----------------
 1 file changed, 10 insertions(+), 17 deletions(-)

diff --git a/src/pspm_sf.m b/src/pspm_sf.m
index 0cad00a31..99ea03a60 100644
--- a/src/pspm_sf.m
+++ b/src/pspm_sf.m
@@ -117,14 +117,7 @@
   warning('ID:number_of_elements_dont_match',...
     'Number of data files and model files does not match.'); return;
 end
-% 2.5 check substhresh --
-if ~isfield(model, 'substhresh')
-  model.substhresh = 2;
-elseif ~isnumeric(model.substhresh)
-  warning('ID:invalid_input', 'Subsession threshold must be numeric.');
-  return;
-end
-% 2.6 check methods --
+% 2.5 check methods --
 if ~isfield(model, 'method')
   model.method = {'dcm'};
 elseif ischar(model.method)
@@ -163,7 +156,7 @@
     end
   end
 end
-% 2.7 Check timing --
+% 2.6 Check timing --
 if strcmpi(model.timeunits, 'whole')
   epochs = repmat({[1 1]}, numel(model.datafile), 1);
 else
@@ -175,19 +168,19 @@
     end
   end
 end
-% 2.8 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.9 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.10 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)
@@ -224,10 +217,10 @@
       data{1}.header.units);
   end
   % 3.5 Get missing epochs --
-  % 3.5.1 Load 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 ---
+  % 3.5.2 sort missing epochs --
     if size(missing{iSn}, 1) > 0
       [~, sortindx] = sort(missing{iSn}(:, 1));
       missing{iSn} = missing{iSn}(sortindx,:);
@@ -242,7 +235,7 @@
   else
     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);
@@ -292,7 +285,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{iSn}(iEpoch, :));
       sf.model{k}(iEpoch).timeunits  = model.timeunits;
@@ -305,7 +298,7 @@
         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
+      % 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

From 0d40d8b14e1caa425f68c1b375a8d423e81a8c85 Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Mon, 3 Jul 2023 01:25:25 +0100
Subject: [PATCH 23/48] Update pspm_load_data.m

---
 src/pspm_load_data.m | 35 ++++++++++-------------------------
 1 file changed, 10 insertions(+), 25 deletions(-)

diff --git a/src/pspm_load_data.m b/src/pspm_load_data.m
index 1ffbc574f..181915721 100644
--- a/src/pspm_load_data.m
+++ b/src/pspm_load_data.m
@@ -98,7 +98,6 @@
     warning('ID:invalid_input', 'Too many inputs specified.');
     return
 end
-
 %% 3 Check fn
 % fn has to be a file or a struct
 switch class(fn)
@@ -153,7 +152,6 @@
     warning('ID:invalid_input', 'fn needs to be an existing file or a struct.');
     return
 end
-
 %% 4 Check channel
 switch class(channel)
   case 'double'
@@ -190,7 +188,6 @@
   otherwise
     warning('ID:invalid_input', 'Unknown channel option.');
 end
-
 %% 5 Check infos
 if isstruct(channel)
   infos = channel.infos;
@@ -217,8 +214,7 @@
   warning('ID:invalid_data_structure', 'Input data does not have sufficient infos');
   return
 end
-
-%% 6 Check data
+%% 6 Load data
 if isstruct(channel)
   data = channel.data;
 else
@@ -232,6 +228,7 @@
     clear loaded_data
   end
 end
+%% 7 Check data
 % initialise error flags
 vflag = zeros(numel(data), 1); % records data structure, valid if 0
 wflag = zeros(numel(data), 1); % records whether data is out of range, valid if 0
@@ -239,7 +236,12 @@
 zflag = zeros(numel(data), 1); % records whether data is empty
 % loop through channels
 for k = 1:numel(data)
-  % check header
+  % 8.1 Convert header channeltype into chantype if there are --
+  if isfield(data{k}.header, 'channeltype')
+    data{k}.header.chantype = data{k}.header.channeltype;
+    data{k}.header = rmfield( data{k}.header, 'channeltype');
+  end
+  % 8.2 Check header --
   if ~isfield(data{k}, 'header')
     vflag(k) = 1;
   else
@@ -259,7 +261,7 @@
       end
     end
   end
-  % check data
+  % 8.3 Check data --
   if vflag(k)==0 && nflag(k)==0 && flag_infos==0
     % required information is available and valid in header and infos
     if ~isfield(data{k}, 'data')
@@ -307,18 +309,6 @@
   warning('ID:missing_data', 'Channel %01.0f is empty.', find(zflag,1));
   % if there is empty data, give a warning but do not suspend
 end
-
-
-%% 7 Autofill information in header
-% some other optional fields which can be autofilled with default values
-% should be added here.
-for k = 1:numel(data)
-  if isfield(data{k}.header, 'channeltype')
-    data{k}.header.chantype = data{k}.header.channeltype;
-    data{k}.header = rmfield( data{k}.header, 'channeltype');
-  end
-end
-
 %% 8 Analyse file structure
 filestruct.numofwavechan = 0;
 filestruct.numofeventchan = 0;
@@ -339,7 +329,6 @@
 elseif numel(filestruct.posofmarker) > 1
   filestruct.posofmarker = filestruct.posofmarker(1); % first marker channel
 end
-
 %% 9 Return channels, or save file
 if isstruct(channel)
   infos = channel.infos;
@@ -391,7 +380,6 @@
 else
   filestruct.posofchannels = [];
 end
-
 sts = 1;
 return
 
@@ -438,7 +426,6 @@
 combined_channels = combined_channels & pupil_channels;
 besteye_channels = besteye_channels & pupil_channels & ~preprocessed_channels;
 % best eye will not select preprocessed eyes
-
 if any(combined_channels)
   flag = combined_channels;
 elseif any(preprocessed_channels) && ~prefer_unprocessed
@@ -450,7 +437,6 @@
   flag = besteye_channels;
 end
 
-
 function flag = get_chans_to_load_for_sps(data, best_eye)
 % 16-06-21 This is a tempory patch for loading sps data, copied from
 % pupil data
@@ -494,7 +480,6 @@
 preprocessed_channels = preprocessed_channels & sps_channels;
 combined_channels = combined_channels & sps_channels;
 besteye_channels = besteye_channels & sps_channels;
-
 if any(combined_channels)
   flag = combined_channels;
 elseif any(preprocessed_channels)
@@ -504,4 +489,4 @@
   end
 else
   flag = besteye_channels;
-end
+end
\ No newline at end of file

From 1e81f0a014c9216c5f60d93d77e4590e220dde87 Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Mon, 3 Jul 2023 01:27:35 +0100
Subject: [PATCH 24/48] Update pspm_load_data.m

---
 src/pspm_load_data.m | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/pspm_load_data.m b/src/pspm_load_data.m
index 181915721..8dc30004d 100644
--- a/src/pspm_load_data.m
+++ b/src/pspm_load_data.m
@@ -229,22 +229,22 @@
   end
 end
 %% 7 Check data
-% initialise error flags
+% 7.1 initialise error flags --
 vflag = zeros(numel(data), 1); % records data structure, valid if 0
 wflag = zeros(numel(data), 1); % records whether data is out of range, valid if 0
 nflag = zeros(numel(data), 1);
 zflag = zeros(numel(data), 1); % records whether data is empty
 % loop through channels
 for k = 1:numel(data)
-  % 8.1 Convert header channeltype into chantype if there are --
-  if isfield(data{k}.header, 'channeltype')
-    data{k}.header.chantype = data{k}.header.channeltype;
-    data{k}.header = rmfield( data{k}.header, 'channeltype');
-  end
-  % 8.2 Check header --
+  % 7.2 Check header --
   if ~isfield(data{k}, 'header')
     vflag(k) = 1;
   else
+    % 7.2.1 Convert header channeltype into chantype if there are --
+    if isfield(data{k}.header, 'channeltype')
+      data{k}.header.chantype = data{k}.header.channeltype;
+      data{k}.header = rmfield(data{k}.header, 'channeltype');
+    end
     if (~isfield(data{k}.header, 'channeltype') && ~isfield(data{k}.header, 'chantype')) || ...
         ~isfield(data{k}.header, 'sr') || ...
         ~isfield(data{k}.header, 'units')
@@ -261,7 +261,7 @@
       end
     end
   end
-  % 8.3 Check data --
+  % 7.3 Check data --
   if vflag(k)==0 && nflag(k)==0 && flag_infos==0
     % required information is available and valid in header and infos
     if ~isfield(data{k}, 'data')

From 0b8e1af6ba335b654683707c79e4262567579262 Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Mon, 3 Jul 2023 01:39:51 +0100
Subject: [PATCH 25/48] Update pspm_load_data.m

---
 src/pspm_load_data.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pspm_load_data.m b/src/pspm_load_data.m
index 8dc30004d..c0d58822e 100644
--- a/src/pspm_load_data.m
+++ b/src/pspm_load_data.m
@@ -245,7 +245,7 @@
       data{k}.header.chantype = data{k}.header.channeltype;
       data{k}.header = rmfield(data{k}.header, 'channeltype');
     end
-    if (~isfield(data{k}.header, 'channeltype') && ~isfield(data{k}.header, 'chantype')) || ...
+    if ~isfield(data{k}.header, 'chantype') || ...
         ~isfield(data{k}.header, 'sr') || ...
         ~isfield(data{k}.header, 'units')
       vflag(k) = 1;

From dea31ee8c36dfc7eb96d127ba3633aac102816e5 Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Mon, 3 Jul 2023 03:05:42 +0100
Subject: [PATCH 26/48] Update pspm_sf_dcm.m

---
 src/pspm_sf_dcm.m | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/pspm_sf_dcm.m b/src/pspm_sf_dcm.m
index 688b59ba8..93b2ab0d2 100644
--- a/src/pspm_sf_dcm.m
+++ b/src/pspm_sf_dcm.m
@@ -24,7 +24,9 @@
 %   │              power, a sensible size is 60 s at 10 Hz)
 %   ├────────.sr:  sampling rate in Hz
 %   │ ▶︎ Optional
-%   └.missing_data: the datafile of missing epochs
+%   └.missing_data: missing epoch data, originally loaded as model.missing
+%                  from pspm_sf, but calculated into .missing_data (created
+%                  in pspm_sf and then transferred to pspm_sf_dcm.
 %
 %   ┌────options:  options structure
 %   ├─.threshold:  threshold for SN detection (default 0.1 mcS)

From ed24fe461c161d11cab4d3e484ad2df328a234c9 Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Mon, 3 Jul 2023 10:35:35 +0100
Subject: [PATCH 27/48] Update pspm_sf_dcm.m

---
 src/pspm_sf_dcm.m | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/pspm_sf_dcm.m b/src/pspm_sf_dcm.m
index 93b2ab0d2..bd6ff284c 100644
--- a/src/pspm_sf_dcm.m
+++ b/src/pspm_sf_dcm.m
@@ -24,7 +24,7 @@
 %   │              power, a sensible size is 60 s at 10 Hz)
 %   ├────────.sr:  sampling rate in Hz
 %   │ ▶︎ Optional
-%   └.missing_data: missing epoch data, originally loaded as model.missing
+%   └.missing_data:missing epoch data, originally loaded as model.missing
 %                  from pspm_sf, but calculated into .missing_data (created
 %                  in pspm_sf and then transferred to pspm_sf_dcm.
 %
@@ -37,7 +37,7 @@
 %   ├.dispsmallwin:
 %   │              display intermediate windows (default 0);
 %   └─.missingthresh:
-%               threshold value for controlling missing epochs (default 2s).
+%                  threshold value for controlling missing epochs (default 2s).
 % ● References
 %   Bach DR, Daunizeau J, Kuelzow N, Friston KJ, & Dolan RJ (2011). Dynamic
 %   causal modelling of spontaneous fluctuations in skin conductance.
@@ -58,7 +58,7 @@
 % 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
+try model.missing_data; catch, warning('Missing epoch data index is not loaded.'); end
 % 2.2 Validate parameters ---
 if ~isnumeric(model.sr) || numel(model.sr) > 1
   errmsg = sprintf('No valid sample rate given.');

From 7072b7216de5d6d930eb28d5db65043a0ced9a71 Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Mon, 3 Jul 2023 11:36:31 +0100
Subject: [PATCH 28/48] update interpolation

---
 src/pspm_guide.m    |  6 ++--
 src/pspm_interp1.m  | 84 +++++++++++++++++++++++++++++++++++++++++++++
 src/pspm_prepdata.m | 60 --------------------------------
 src/pspm_sf_dcm.m   |  6 ++--
 4 files changed, 91 insertions(+), 65 deletions(-)
 create mode 100644 src/pspm_interp1.m

diff --git a/src/pspm_guide.m b/src/pspm_guide.m
index da1a514a7..8c95bba8f 100644
--- a/src/pspm_guide.m
+++ b/src/pspm_guide.m
@@ -1,7 +1,7 @@
 function varargout = pspm_guide(varargin)
-% �? Description
+% ● Description
 %   pspm_guide handles the main GUI for PsPM
-% �? Developer's Guide
+% ● Developer's Guide
 %   * Template
 %     function varargout = FunctionName(hObject, eventdata, handles, varargin)
 %       varargout  cell array for returning output args (see VARARGOUT);
@@ -15,7 +15,7 @@
 %       contents{get(hObject,'Value')} returns selected item from the list
 %   * Hint
 %       popupmenu controls usually have a white background on Windows.
-% �? Copyright
+% ● Copyright
 %   Introduced in PsPM 1.0 and terminated in PsPM 6.1.
 %   Written by Dominik R Bach (Wellcome Centre for Human Neuroimaging) in 2008-2021
 %   Lastly updated in PsPM 6.1 by Teddy Chao (UCL) in 2022
diff --git a/src/pspm_interp1.m b/src/pspm_interp1.m
new file mode 100644
index 000000000..f4dd7b467
--- /dev/null
+++ b/src/pspm_interp1.m
@@ -0,0 +1,84 @@
+function Y = pspm_interp1(varargin)
+% ● Description
+%   pspm_interp1 is a shared PsPM function for interpolating data with NaNs
+%   based on the reference of missing epochs and first order interpolation.
+% ● Format
+%   Y = pspm_interp1(varargin)
+% ● Arguments
+%   X:              data that contains NaNs to be interpolated
+%   index_missing:  index of missing epochs with the same size of X in binary
+%                   values. 1 if NaNs, 0 if non-NaNs.
+%   Y:              processed data
+% ● History
+%   Introduced In PsPM 6.1
+%   Written in 2023 by Teddy Chao (UCL)
+
+switch nargin
+  case 1
+    X = varargin{1};
+    index_missing = zeros(size(X));
+  case 2
+    X = varargin{1};
+    index_missing = varargin{2};
+  otherwise
+    warning('ID:invalid_input','pspm_interp1 accepts up to two arguments');
+end
+% find nan head and tail
+X_nan_head = 0;
+X_nan_tail = 0;
+X_nan_head_range = [];
+X_nan_tail_range = [];
+X_nan_head_interp = [];
+X_nan_tail_interp = [];
+index_non_nan_full = 1:length(X);
+index_non_nan_full = index_non_nan_full(~isnan(X));
+if index_non_nan_full(1) > 1
+  X_nan_head = 1;
+  X_nan_head_range = 1:(index_non_nan_full(1)-1);
+end
+if index_non_nan_full(end) < length(X)
+  X_nan_tail = 1;
+  X_nan_tail_range = (index_non_nan_full(end)+1):length(X);
+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);
+else
+  X_body_interp = X_body;
+end
+% interpolate head
+if X_nan_head
+  X_nan_head_interp = interp1(...
+    (1:length(X_body_interp))+length(X_nan_head_range),...
+    X_body_interp,...
+    X_nan_head_range,'linear','extrap');
+end
+% interpolate tail
+if X_nan_tail
+  X_nan_tail_interp = interp1(...
+    (1:length(X_body_interp))+length(X_nan_head_range),...
+    X_body_interp,...
+    X_nan_tail_range,'linear','extrap');
+end
+if iscolumn(X_body_interp)
+  if ~iscolumn(X_nan_head_interp)
+    X_nan_head_interp = transpose(X_nan_head_interp);
+  end
+  if ~iscolumn(X_nan_tail_interp)
+    X_nan_tail_interp = transpose(X_nan_tail_interp);
+  end
+  Y = [X_nan_head_interp; X_body_interp; X_nan_tail_interp];
+else
+  if iscolumn(X_nan_head_interp)
+    X_nan_head_interp = transpose(X_nan_head_interp);
+  end
+  if iscolumn(X_nan_tail_interp)
+    X_nan_tail_interp = transpose(X_nan_tail_interp);
+  end
+  Y = [X_nan_head_interp, X_body_interp, X_nan_tail_interp];
+end
+return
diff --git a/src/pspm_prepdata.m b/src/pspm_prepdata.m
index a2838a471..d9b2a4f0b 100644
--- a/src/pspm_prepdata.m
+++ b/src/pspm_prepdata.m
@@ -197,63 +197,3 @@
     varargout{3} = newsr;
 end
 return
-function data_interp = pspm_interp1(data)
-% find nan head and tail
-data_nan_head = 0;
-data_nan_tail = 0;
-data_nan_head_range = [];
-data_nan_tail_range = [];
-data_nan_head_interp = [];
-data_nan_tail_interp = [];
-index_non_nan_full = 1:length(data);
-index_non_nan_full = index_non_nan_full(~isnan(data));
-if index_non_nan_full(1) > 1
-  data_nan_head = 1;
-  data_nan_head_range = 1:(index_non_nan_full(1)-1);
-end
-if index_non_nan_full(end) < length(data)
-  data_nan_tail = 1;
-  data_nan_tail_range = (index_non_nan_full(end)+1):length(data);
-end
-data_body = data(index_non_nan_full(1):index_non_nan_full(end));
-% processing body
-index = 1:length(data_body);
-index_nan = index(isnan(data_body));
-index_non_nan = index(~isnan(data_body));
-if ~isempty(index_nan)
-  data_body_interp = interp1(index_non_nan,data_body(index_non_nan),index_nan);
-else
-  data_body_interp = data_body;
-end
-% interpolate head
-if data_nan_head
-  data_nan_head_interp = interp1(...
-    (1:length(data_body_interp))+length(data_nan_head_range),...
-    data_body_interp,...
-    data_nan_head_range,'linear','extrap');
-end
-% interpolate tail
-if data_nan_tail
-  data_nan_tail_interp = interp1(...
-    (1:length(data_body_interp))+length(data_nan_head_range),...
-    data_body_interp,...
-    data_nan_tail_range,'linear','extrap');
-end
-if iscolumn(data_body_interp)
-  if ~iscolumn(data_nan_head_interp)
-    data_nan_head_interp = transpose(data_nan_head_interp);
-  end
-  if ~iscolumn(data_nan_tail_interp)
-    data_nan_tail_interp = transpose(data_nan_tail_interp);
-  end
-  data_interp = [data_nan_head_interp; data_body_interp; data_nan_tail_interp];
-else
-  if iscolumn(data_nan_head_interp)
-    data_nan_head_interp = transpose(data_nan_head_interp);
-  end
-  if iscolumn(data_nan_tail_interp)
-    data_nan_tail_interp = transpose(data_nan_tail_interp);
-  end
-  data_interp = [data_nan_head_interp, data_body_interp, data_nan_tail_interp];
-end
-return
diff --git a/src/pspm_sf_dcm.m b/src/pspm_sf_dcm.m
index bd6ff284c..17ad34b1b 100644
--- a/src/pspm_sf_dcm.m
+++ b/src/pspm_sf_dcm.m
@@ -162,9 +162,11 @@
   flag_missing_too_long = 1;
 end
 options.isYout = ymissing(:)';
+% 4.6 interpolate data body to fill NaNs
+y_interpolated = pspm_interp1(y, ymissing);
 %% 5 Extract parameters
 if ~flag_missing_too_long
-  [posterior, output] = VBA_NLStateSpaceModel(y(:)',u,f_fname,g_fname,dim,options);
+  [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
@@ -188,7 +190,7 @@
   out.model.posterior = posterior;
   out.model.output    = output;
   out.model.u         = u;
-  out.model.y         = y(:)';
+  out.model.y         = y_interpolated(:)';
   out.time            = toc(tstart);
 else
   out.t               = NaN;

From 8d5af99cf8cfd7cfbcac423d5954bb6370a8b0a6 Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Mon, 3 Jul 2023 11:41:06 +0100
Subject: [PATCH 29/48] Update pspm_interp1.m

---
 src/pspm_interp1.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pspm_interp1.m b/src/pspm_interp1.m
index f4dd7b467..bf74583d3 100644
--- a/src/pspm_interp1.m
+++ b/src/pspm_interp1.m
@@ -10,7 +10,7 @@
 %                   values. 1 if NaNs, 0 if non-NaNs.
 %   Y:              processed data
 % ● History
-%   Introduced In PsPM 6.1
+%   Introduced in PsPM 6.1
 %   Written in 2023 by Teddy Chao (UCL)
 
 switch nargin

From 989f791410e73ce3f4ce7e8b50d287bd258f155d Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Mon, 3 Jul 2023 13:18:43 +0100
Subject: [PATCH 30/48] update sf options

---
 src/pspm_options.m | 2 +-
 src/pspm_sf.m      | 4 +++-
 src/pspm_sf_dcm.m  | 9 +++++----
 3 files changed, 9 insertions(+), 6 deletions(-)

diff --git a/src/pspm_options.m b/src/pspm_options.m
index 09ee2086d..ce579b88a 100644
--- a/src/pspm_options.m
+++ b/src/pspm_options.m
@@ -413,6 +413,7 @@
     options = autofill(options,'marker_chan_num',         1,          '*Int*Char'       );
     options = autofill(options,'overwrite',               1,          0                 );
     options = autofill(options,'threshold',               0.1,        '>', 0            );
+    options = autofill(options,'missingthresh',           2,          '>', 0            );
     options = autofill(options,'theta',                   [0.923581, ...
                                                           3.921034, ...
                                                           2.159389, ...
@@ -424,7 +425,6 @@
     options = autofill(options,'dispwin',                 1,          0                 );
     options = autofill(options,'dispsmallwin',            0,          1                 );
     options = autofill(options,'fresp',                   0.5,        '>', 0            );
-    options = autofill(options,'missingthresh',           2,          '>', 0            );
     options = autofill(options,'threshold',               0.1,        '>', 0            );
     options = autofill(options,'theta',                   [0.923581, ...
                                                           3.921034, ...
diff --git a/src/pspm_sf.m b/src/pspm_sf.m
index 384b6d4cf..7c590df8c 100644
--- a/src/pspm_sf.m
+++ b/src/pspm_sf.m
@@ -45,7 +45,9 @@
 %   │                   (default: read from pspm_sf_theta)
 %   ├──────────.fresp:  frequency of responses to model (default 0.5 Hz)
 %   ├────────.dispwin:  display progress window (default 1)
-%   └───.dispsmallwin:  display intermediate windows (default 0);
+%   ├───.dispsmallwin:  display intermediate windows (default 0);
+%   └─.missingthresh:
+%                       threshold value for controlling missing epochs (default 2s).
 % ● References
 %   1.[DCM for SF]
 %     Bach DR, Daunizeau J, Kuelzow N, Friston KJ, Dolan RJ (2010). Dynamic
diff --git a/src/pspm_sf_dcm.m b/src/pspm_sf_dcm.m
index 17ad34b1b..46e3def79 100644
--- a/src/pspm_sf_dcm.m
+++ b/src/pspm_sf_dcm.m
@@ -36,8 +36,9 @@
 %   ├───.dispwin:  display progress window (default 1)
 %   ├.dispsmallwin:
 %   │              display intermediate windows (default 0);
-%   └─.missingthresh:
-%                  threshold value for controlling missing epochs (default 2s).
+%   └.missingthresh:
+%                  threshold value for controlling missing epochs (default 2s),
+%                  which is originally inherited from SF
 % ● References
 %   Bach DR, Daunizeau J, Kuelzow N, Friston KJ, & Dolan RJ (2011). Dynamic
 %   causal modelling of spontaneous fluctuations in skin conductance.
@@ -83,8 +84,8 @@
 if options.invalid
   return
 end
-options.DisplayWin = options.dispwin;
-options.GnFigs = options.dispsmallwin;
+% options.DisplayWin = options.dispwin;
+% options.GnFigs = options.dispsmallwin;
 fresp = options.fresp;
 threshold = options.threshold;
 try

From 4dcf6a8430fa5c5ff2c189764bb58fe11afa896c Mon Sep 17 00:00:00 2001
From: "teddy.chao" <dadi.zhao@ucl.ac.uk>
Date: Mon, 3 Jul 2023 14:02:17 +0100
Subject: [PATCH 31/48] update sf help text

---
 src/pspm_sf.m     | 34 +++++++++++++++++++---------------
 src/pspm_sf_dcm.m | 25 +++++++++++++++----------
 2 files changed, 34 insertions(+), 25 deletions(-)

diff --git a/src/pspm_sf.m b/src/pspm_sf.m
index 7c590df8c..37ab7b9fc 100644
--- a/src/pspm_sf.m
+++ b/src/pspm_sf.m
@@ -25,29 +25,33 @@
 %   │                   as 'dcm'
 %   │                   [cell_array] a cell array of methods mentioned above.
 %   ├─────────.filter:  filter settings; modality specific default
-%   ├────────.missing:  [string/cell_array]
+%   ├────────.missing:  [string/cell_array] [default: no missing values]
 %   │                   Allows to specify missing (e.g. artefact) epochs in the
 %   │                   data file. See pspm_get_timing for epoch definition; specify
 %   │                   a cell array for multiple input files. This must always be
 %   │                   specified in SECONDS.
-%   │                   Default: no missing values
-%   └────────.channel:  channel number; default: first SCR channel
+%   └────────.channel:  [integer] [default: first SCR channel]
+%                       channel number.
 %   ┌─────────options
-%   ├──────.overwrite:  [logical] (0 or 1)
+%   ├──────.overwrite:  [logical] [default: determined by pspm_overwrite]
 %   │                   Define whether to overwrite existing output files or not.
-%   │                   Default value: determined by pspm_overwrite.
-%   ├.marker_chan_num:  marker channel number
+%   ├.marker_chan_num:  [integer] 
+%   │                   marker channel number
 %   │                   if undefined or 0, first marker channel is used.
-%   │ additional options for individual methods:
+%   │ * Additional options for individual methods:
 %   │ dcm related options
-%   ├──────.threshold:  threshold for SN detection (default 0.1 mcS)
-%   ├──────────.theta:  a (1 x 5) vector of theta values for f_SF
-%   │                   (default: read from pspm_sf_theta)
-%   ├──────────.fresp:  frequency of responses to model (default 0.5 Hz)
-%   ├────────.dispwin:  display progress window (default 1)
-%   ├───.dispsmallwin:  display intermediate windows (default 0);
-%   └─.missingthresh:
-%                       threshold value for controlling missing epochs (default 2s).
+%   ├──────.threshold:  [numeric] [default: 0.1] [unit: mcS]
+%   │                   threshold for SN detection (default 0.1 mcS)
+%   ├──────────.theta:  [vector] [default: read from pspm_sf_theta]
+%   │                   A (1 x 5) vector of theta values for f_SF.
+%   ├──────────.fresp:  [numeric] [unit: Hz] [default: 0.5]
+%   │                   frequency of responses to model.
+%   ├────────.dispwin:  [logical] [default: 1]
+%   │                   display progress window.
+%   ├───.dispsmallwin:  [logical] [default: 0]
+%   │                   display intermediate windows.
+%   └──.missingthresh:  [numeric] [default: 2] [unit: second]
+%                       threshold value for controlling missing epochs.
 % ● References
 %   1.[DCM for SF]
 %     Bach DR, Daunizeau J, Kuelzow N, Friston KJ, Dolan RJ (2010). Dynamic
diff --git a/src/pspm_sf_dcm.m b/src/pspm_sf_dcm.m
index 46e3def79..1b73ceb72 100644
--- a/src/pspm_sf_dcm.m
+++ b/src/pspm_sf_dcm.m
@@ -22,22 +22,27 @@
 %   │ ▶︎ Mandatory
 %   ├───────.scr:  skin conductance epoch (maximum size depends on computing
 %   │              power, a sensible size is 60 s at 10 Hz)
-%   ├────────.sr:  sampling rate in Hz
+%   ├────────.sr:  [numeric] [unit: Hz]
+%   │              sampling rate.
 %   │ ▶︎ Optional
 %   └.missing_data:missing epoch data, originally loaded as model.missing
 %                  from pspm_sf, but calculated into .missing_data (created
 %                  in pspm_sf and then transferred to pspm_sf_dcm.
 %
-%   ┌────options:  options structure
-%   ├─.threshold:  threshold for SN detection (default 0.1 mcS)
-%   ├─────.theta:  a (1 x 5) vector of theta values for f_SF
-%   │              (default: read from pspm_sf_theta)
-%   ├─────.fresp:  frequency of responses to model (default 0.5 Hz)
-%   ├───.dispwin:  display progress window (default 1)
-%   ├.dispsmallwin:
-%   │              display intermediate windows (default 0);
+%   ┌────options
+%   ├─.threshold:  [numeric] [default: 0.1] [unit: mcS]
+%   │              threshold for SN detection (default 0.1 mcS)
+%   ├─────.theta:  [vector] [default: read from pspm_sf_theta]
+%   │               a (1 x 5) vector of theta values for f_SF
+%   ├─────.fresp:  [numeric] [unit: Hz] [default: 0.5]
+%   │              frequency of responses to model
+%   ├───.dispwin:  [logical] [default: 1]
+%   │              display progress window.
+%   ├.dispsmallwin:[logical] [default: 0]
+%   │              display intermediate windows.
 %   └.missingthresh:
-%                  threshold value for controlling missing epochs (default 2s),
+%                  [numeric] [default: 2] [unit: second]
+%                  threshold value for controlling missing epochs,
 %                  which is originally inherited from SF
 % ● References
 %   Bach DR, Daunizeau J, Kuelzow N, Friston KJ, & Dolan RJ (2011). Dynamic

From a12dafa991a203bc25dfadc9a25b0aafb701f7dd Mon Sep 17 00:00:00 2001
From: "teddy.chao" <dadi.zhao@ucl.ac.uk>
Date: Mon, 3 Jul 2023 15:23:21 +0100
Subject: [PATCH 32/48] Update the GUI

---
 src/pspm_cfg/pspm_cfg_run_sf.m |  4 +++-
 src/pspm_cfg/pspm_cfg_sf.m     | 16 ++++++++++++----
 src/pspm_sf_dcm.m              |  1 -
 3 files changed, 15 insertions(+), 6 deletions(-)

diff --git a/src/pspm_cfg/pspm_cfg_run_sf.m b/src/pspm_cfg/pspm_cfg_run_sf.m
index 2da224e68..1a346d0d9 100644
--- a/src/pspm_cfg/pspm_cfg_run_sf.m
+++ b/src/pspm_cfg/pspm_cfg_run_sf.m
@@ -77,7 +77,9 @@
 if ~isempty(job.fresp)
     options.fresp = job.fresp;
 end
-
+if ~isempty(job.missing)
+  options.missing = job.missing;
+end
 options.dispwin = job.dispwin;
 options.dispsmallwin = job.dispsmallwin;
 
diff --git a/src/pspm_cfg/pspm_cfg_sf.m b/src/pspm_cfg/pspm_cfg_sf.m
index fa6411f93..cdc9c0532 100644
--- a/src/pspm_cfg/pspm_cfg_sf.m
+++ b/src/pspm_cfg/pspm_cfg_sf.m
@@ -143,7 +143,7 @@
 direction.labels  = {'Unidirectional', 'Bidirectional'};
 direction.values  = {'uni', 'bi'};
 direction.help    = {['A unidirectional filter is applied twice in the forward direction. ' ...
-    'A �bidirectional� filter is applied once in the forward direction and once in the ' ...
+    'A "bidirectional" filter is applied once in the forward direction and once in the ' ...
     'backward direction to correct the temporal shift due to filtering in forward direction.']};
 
 filter_edit        = cfg_branch;
@@ -228,8 +228,8 @@
 timeunits.tag     = 'timeunits';
 timeunits.values  = {seconds, samples, markers, whole};
 timeunits.help    = {['Indicate the time units on which the specification of the conditions will be based. ' ...
-    'Time units can be specified in �seconds�, number of �markers�, or number of data �samples� . Time units ' ...
-    'refer to the beginning of the data file and not to the beginning of the original recordings e. g. if ' ...
+    'Time units can be specified in "seconds", number of "markers", or number of data "samples". Time units ' ...
+    'refer to the beginning of the data file and not to the beginning of the original recordings e.g. if ' ...
     'data were trimmed.']};
 
 %% Channel nr
@@ -292,6 +292,14 @@
 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',...
+                  'data must always be specified in seconds.'};
+
 % Show figures
 dispwin         = cfg_menu;
 dispwin.name    = 'Display Progress Window';
@@ -314,7 +322,7 @@
 sf      = cfg_exbranch;
 sf.name = 'SF';
 sf.tag  = 'sf';
-sf.val  = {datafile, modelfile, outdir, method, timeunits, filter, chan, overwrite, threshold, theta, fresp, dispwin, dispsmallwin};
+sf.val  = {datafile, modelfile, outdir, method, timeunits, filter, chan, overwrite, threshold, missing, theta, fresp, dispwin, dispsmallwin};
 sf.prog = @pspm_cfg_run_sf;
 sf.vout = @pspm_cfg_vout_sf;
 sf.help = {['This suite of models is designed for analysing spontaneous fluctuations (SF) in skin conductance ' ...
diff --git a/src/pspm_sf_dcm.m b/src/pspm_sf_dcm.m
index 1b73ceb72..0f2b5c510 100644
--- a/src/pspm_sf_dcm.m
+++ b/src/pspm_sf_dcm.m
@@ -64,7 +64,6 @@
 % 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 epoch data index is not loaded.'); end
 % 2.2 Validate parameters ---
 if ~isnumeric(model.sr) || numel(model.sr) > 1
   errmsg = sprintf('No valid sample rate given.');

From 6920e19f0dcbafa74a7321518fa352b9f4fc1c70 Mon Sep 17 00:00:00 2001
From: "teddy.chao" <dadi.zhao@ucl.ac.uk>
Date: Mon, 3 Jul 2023 15:28:26 +0100
Subject: [PATCH 33/48] Update pspm_cfg_sf.m

---
 src/pspm_cfg/pspm_cfg_sf.m | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/pspm_cfg/pspm_cfg_sf.m b/src/pspm_cfg/pspm_cfg_sf.m
index cdc9c0532..ab22fc507 100644
--- a/src/pspm_cfg/pspm_cfg_sf.m
+++ b/src/pspm_cfg/pspm_cfg_sf.m
@@ -297,8 +297,8 @@
 missing.tag     = 'missingfile';
 missing.num     = [1 1];
 missing.filter  = '.*\.(mat|MAT)$';
-missing.help    = {'Missing (e.g. artefact) epochs in the data file, where',...
-                  'data must always be specified in seconds.'};
+missing.help    = {['Missing (e.g. artefact) epochs in the data file, where ',...
+                  'data must always be specified in seconds.']};
 
 % Show figures
 dispwin         = cfg_menu;

From b29e6cc2a9402fede73949da7e65aaa95fa5cf30 Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Fri, 7 Jul 2023 12:58:13 +0100
Subject: [PATCH 34/48] Update pspm_interp1.m

---
 src/pspm_interp1.m | 29 ++++++++++++++++++++++++++++-
 1 file changed, 28 insertions(+), 1 deletion(-)

diff --git a/src/pspm_interp1.m b/src/pspm_interp1.m
index bf74583d3..f0c8b827e 100644
--- a/src/pspm_interp1.m
+++ b/src/pspm_interp1.m
@@ -13,6 +13,7 @@
 %   Introduced in PsPM 6.1
 %   Written in 2023 by Teddy Chao (UCL)
 
+%% 1 Load inputs
 switch nargin
   case 1
     X = varargin{1};
@@ -23,7 +24,33 @@
   otherwise
     warning('ID:invalid_input','pspm_interp1 accepts up to two arguments');
 end
-% find nan head and tail
+%% 2 Check inputs
+switch sum(~isnan(a))
+  case 0
+    % if there are no non-nans, do not process any interpolation, give a
+    % warning and return
+    warning('ID:invalid_input',...
+      'Input data contains only NaNs thus cannot be interpolated.')
+    Y = X;
+    return
+  case 1
+    % if there are only 1 non-nan, do not process any interpolation,
+    % give a warning and explain the reason
+    warning('ID:invalid_input',...
+      'Input data contains only 1 non-NaN thus cannot be interpolated.')
+    Y = X;
+    return
+  otherwise
+    % if there are less than 10^ non-nan, still perform interpolation,
+    % however give a warning and explain the reason
+    non_nan_percentage = sum(~isnan(a))/length(a);
+    if non_nan_percentage<0.1
+      warning('ID:invalid_input',...
+      'Input data contains less than 10% non-NaN. Interpolation can ',... 
+      'still be performed but results could be inaccurate.')
+    end
+end
+%% 3 find nan head and tail
 X_nan_head = 0;
 X_nan_tail = 0;
 X_nan_head_range = [];

From 29f951f2cc1fdd94cc0a2a171b7bb7d357967504 Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Fri, 7 Jul 2023 13:17:52 +0100
Subject: [PATCH 35/48] Update pspm_interp1.m

---
 src/pspm_interp1.m | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/pspm_interp1.m b/src/pspm_interp1.m
index f0c8b827e..86b1e5aae 100644
--- a/src/pspm_interp1.m
+++ b/src/pspm_interp1.m
@@ -25,7 +25,7 @@
     warning('ID:invalid_input','pspm_interp1 accepts up to two arguments');
 end
 %% 2 Check inputs
-switch sum(~isnan(a))
+switch sum(~isnan(X))
   case 0
     % if there are no non-nans, do not process any interpolation, give a
     % warning and return
@@ -43,7 +43,7 @@
   otherwise
     % if there are less than 10^ non-nan, still perform interpolation,
     % however give a warning and explain the reason
-    non_nan_percentage = sum(~isnan(a))/length(a);
+    non_nan_percentage = sum(~isnan(X))/length(X);
     if non_nan_percentage<0.1
       warning('ID:invalid_input',...
       'Input data contains less than 10% non-NaN. Interpolation can ',... 

From e4986cff1729b6c056df24f2ebd7a44127104859 Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Mon, 10 Jul 2023 13:18:59 +0100
Subject: [PATCH 36/48] Update pspm_resp_pp.m

---
 src/pspm_resp_pp.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pspm_resp_pp.m b/src/pspm_resp_pp.m
index ed3f954cc..96fbb561d 100644
--- a/src/pspm_resp_pp.m
+++ b/src/pspm_resp_pp.m
@@ -85,7 +85,7 @@
 filt.hporder   = 1;
 filt.direction = 'bi';
 filt.down      = 10;
-[sts, newresp, newsr] = pspm_prepdata(resp - mean(resp), filt);
+[sts, newresp, newsr] = pspm_prepdata(resp - mean(resp,"omitnan"), filt);
 % Median filter
 newresp = medfilt1(newresp, ceil(newsr) + 1);
 %% detect breathing cycles

From a57a575af49dbe77785d41ba225d11571c60ec5c Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Thu, 13 Jul 2023 04:16:28 +0100
Subject: [PATCH 37/48] Update pspm_cfg_sf.m

---
 src/pspm_cfg/pspm_cfg_sf.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pspm_cfg/pspm_cfg_sf.m b/src/pspm_cfg/pspm_cfg_sf.m
index ab22fc507..7a890f236 100644
--- a/src/pspm_cfg/pspm_cfg_sf.m
+++ b/src/pspm_cfg/pspm_cfg_sf.m
@@ -195,7 +195,7 @@
 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 ' ...
-    'specified the time units as �markers�.']};
+    'specified the time units as "markers".']};
 
 
 %% Timeunits

From 52d3f99c42f98d67931edaf858a25c8587f62c0f Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Mon, 17 Jul 2023 13:41:04 +0100
Subject: [PATCH 38/48] update sf series functions

---
 src/pspm_cfg/pspm_cfg_run_sf.m |  4 +++-
 src/pspm_cfg/pspm_cfg_sf.m     | 33 +++++++++++++++++++++++++--------
 src/pspm_options.m             |  2 +-
 src/pspm_overwrite.m           |  2 ++
 src/pspm_sf.m                  |  8 ++++----
 src/pspm_sf_mp.m               | 13 +++++++++----
 6 files changed, 44 insertions(+), 18 deletions(-)

diff --git a/src/pspm_cfg/pspm_cfg_run_sf.m b/src/pspm_cfg/pspm_cfg_run_sf.m
index 1a346d0d9..9638d4f9c 100644
--- a/src/pspm_cfg/pspm_cfg_run_sf.m
+++ b/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;
diff --git a/src/pspm_cfg/pspm_cfg_sf.m b/src/pspm_cfg/pspm_cfg_sf.m
index ab22fc507..f6b89ddf6 100644
--- a/src/pspm_cfg/pspm_cfg_sf.m
+++ b/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,11 +193,11 @@
 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 ' ...
-    'specified the time units as �markers�.']};
+    'specified the time units as "markers".']};
 
 
 %% Timeunits
@@ -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';
diff --git a/src/pspm_options.m b/src/pspm_options.m
index ce579b88a..585238b3d 100644
--- a/src/pspm_options.m
+++ b/src/pspm_options.m
@@ -579,7 +579,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')
diff --git a/src/pspm_overwrite.m b/src/pspm_overwrite.m
index f6c446550..f8aa1e642 100644
--- a/src/pspm_overwrite.m
+++ b/src/pspm_overwrite.m
@@ -96,5 +96,7 @@
   case 2
     varargout{1} = sts;
     varargout{2} = overwrite_final;
+  otherwise
+    varargout{1} = overwrite_final;
 end
 return
diff --git a/src/pspm_sf.m b/src/pspm_sf.m
index 37ab7b9fc..e50c94a90 100644
--- a/src/pspm_sf.m
+++ b/src/pspm_sf.m
@@ -201,7 +201,7 @@
   % 3.1 User output
   fprintf('SF analysis: %s ...', model.datafile{iFile});
   % 3.2 Check whether model file exists
-  if ~pspm_overwrite(model.modelfile, options)
+  if ~pspm_overwrite(model.modelfile, options.overwrite)
     return
   end
   % 3.3 get and filter data
@@ -243,7 +243,7 @@
   % 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);
+      [nsts, ~, ndata] = pspm_load_data(model.datafile{iFile}, options.marker_chan_num);
       if nsts == -1
         warning('ID:invalid_input', 'Could not load data');
         return;
@@ -253,14 +253,14 @@
           ['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');
+        [nsts, ~, ~] = pspm_load_data(model.datafile{iFile}, 'marker');
         if nsts == -1
           warning('ID:invalid_input', 'Could not load data');
           return;
         end
       end
     else
-      [nsts, ~, ~] = pspm_load_data(model.datafile, 'marker');
+      [nsts, ~, ~] = pspm_load_data(model.datafile{iFile}, 'marker');
       if nsts == -1
         warning('ID:invalid_input', 'Could not load data');
         return;
diff --git a/src/pspm_sf_mp.m b/src/pspm_sf_mp.m
index 944ae5379..a7544c7bb 100644
--- a/src/pspm_sf_mp.m
+++ b/src/pspm_sf_mp.m
@@ -1,4 +1,4 @@
-function varargout = pspm_sf_mp(scr, sr, options)
+function varargout = pspm_sf_mp(model, options)
 % ● Description
 %   pspm_sf_mp does the inversion of a DCM for SF of the skin conductance, using
 %   a matching pursuit algorithm, and f_SF for the forward model
@@ -47,9 +47,12 @@
 sts = -1;
 tstart = tic;
 
+sr = model.sr;
+scr = model.scr;
+
 % check input arguments
 % ------------------------------------------------------------------------
-if nargin < 2 || ~isnumeric(sr) || numel(sr) > 1
+if ~isnumeric(sr) || numel(sr) > 1
   errmsg = sprintf('No valid sample rate given.');
 elseif (sr < 1) || (sr > 1e5)
   errmsg = sprintf('Sample rate out of range.');
@@ -73,6 +76,7 @@
 S.dt = 1/sr;                    % sampling interval of the data
 S.n = numel(scr);               % n: number of samples in the data segment
 S.sfduration = 30;              % duration of the modelled SF
+S.fresp = options.fresp;
 S.sftail = 10;                  % model tail of SF in previous seconds
 S.ntail = S.sftail/S.dt;        % number of samples in SF tail to model
 S.nsf = S.sfduration/S.dt;      % number of samples in modelled SF
@@ -82,6 +86,7 @@
 S.theta = pspm_sf_theta;         % get SF CRF
 S.maxres = 0.001 * S.n;         % residual threshold per sample
 S.options = options;
+S.threshold = options.threshold;
 
 % generate over-complete dictionary D.D
 % -------------------------------------------------------------------------
@@ -252,9 +257,9 @@
 sts = 1;
 switch nargout
   case 1
-    varargout{1} = mp;
+    varargout{1} = out;
   case 2
     varargout{1} = sts;
-    varargout{2} = mp;
+    varargout{2} = out;
 end
 return

From eb1e6e2376df330891ea970d4397a8f6978eaa25 Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Mon, 17 Jul 2023 14:58:27 +0100
Subject: [PATCH 39/48] update sf marker definition

---
 src/pspm_sf.m     | 2 +-
 src/pspm_sf_auc.m | 6 +++++-
 src/pspm_sf_mp.m  | 4 +++-
 src/pspm_sf_scl.m | 4 +++-
 4 files changed, 12 insertions(+), 4 deletions(-)

diff --git a/src/pspm_sf.m b/src/pspm_sf.m
index e50c94a90..8c855b67e 100644
--- a/src/pspm_sf.m
+++ b/src/pspm_sf.m
@@ -279,7 +279,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
diff --git a/src/pspm_sf_auc.m b/src/pspm_sf_auc.m
index bc2b2a0de..96bdc378a 100644
--- a/src/pspm_sf_auc.m
+++ b/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;
diff --git a/src/pspm_sf_mp.m b/src/pspm_sf_mp.m
index a7544c7bb..8e856b4de 100644
--- a/src/pspm_sf_mp.m
+++ b/src/pspm_sf_mp.m
@@ -47,8 +47,10 @@
 sts = -1;
 tstart = tic;
 
-sr = model.sr;
+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;
 
 % check input arguments
 % ------------------------------------------------------------------------
diff --git a/src/pspm_sf_scl.m b/src/pspm_sf_scl.m
index e82c7c5d4..f8e1d9ee7 100644
--- a/src/pspm_sf_scl.m
+++ b/src/pspm_sf_scl.m
@@ -1,4 +1,4 @@
-function varargout = pspm_sf_scl(scr, sr, options)
+function varargout = pspm_sf_scl(model, options)
 % ● Description
 %   pspm_sf_scl returns the mean skin conductance level for an epoch
 % ● Format
@@ -25,6 +25,8 @@
 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
 scl = mean(scr);
 sts = 1;
 switch nargout

From e35e467522190f9d9370588d9ed9cebf6e782fbc Mon Sep 17 00:00:00 2001
From: "teddy.chao" <dadi.zhao@ucl.ac.uk>
Date: Mon, 24 Jul 2023 02:26:48 +0100
Subject: [PATCH 40/48] Update pspm_sf.m

---
 src/pspm_sf.m | 25 ++++---------------------
 1 file changed, 4 insertions(+), 21 deletions(-)

diff --git a/src/pspm_sf.m b/src/pspm_sf.m
index 8c855b67e..142cb9beb 100644
--- a/src/pspm_sf.m
+++ b/src/pspm_sf.m
@@ -244,29 +244,12 @@
   if any(strcmp(model.timeunits, {'marker', 'markers'}))
     if options.marker_chan_num
       [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
-      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{iFile}, 'marker');
-        if nsts == -1
-          warning('ID:invalid_input', 'Could not load data');
-          return;
-        end
-      end
+      if nsts == -1; warning('ID:invalid_input', 'Could not load data'); return; end
     else
-      [nsts, ~, ~] = pspm_load_data(model.datafile{iFile}, '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

From 9b538c0a0503b5d3efb7ef1b283bc3dc178b4d64 Mon Sep 17 00:00:00 2001
From: "teddy.chao" <dadi.zhao@ucl.ac.uk>
Date: Mon, 24 Jul 2023 02:41:40 +0100
Subject: [PATCH 41/48] Update pspm_sf_mp.m

---
 src/pspm_sf_mp.m | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/pspm_sf_mp.m b/src/pspm_sf_mp.m
index 8e856b4de..6a9f8f0b2 100644
--- a/src/pspm_sf_mp.m
+++ b/src/pspm_sf_mp.m
@@ -164,12 +164,15 @@
 S.nosf = 0;
 k = 1;
 S.Yres = y;
-ind = [];
 % initialise diagnostics
 S.diagnostics.neg = 0; % stop because of zero or negative amplitude
 S.diagnostics.num = 0; % number of iterations
 S.diagnostics.error = NaN; % error
 
+a = [];
+asf = [];
+ind = [];
+
 % run algorithm ---
 while S.cont
   % remove used atoms from dictionary

From 0ef3fdc4cb730d2681cf2fb1eeb5697024609f5f Mon Sep 17 00:00:00 2001
From: "teddy.chao" <dadi.zhao@ucl.ac.uk>
Date: Mon, 24 Jul 2023 14:56:59 +0100
Subject: [PATCH 42/48] Update pspm_sf.m

---
 src/pspm_sf.m | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/pspm_sf.m b/src/pspm_sf.m
index 06df85c9d..2c472ebed 100644
--- a/src/pspm_sf.m
+++ b/src/pspm_sf.m
@@ -244,10 +244,10 @@
   % 3.6 Get marker data --
   if any(strcmp(model.timeunits, {'marker', 'markers'}))
     if options.marker_chan_num
-      [nsts, ~, ndata] = pspm_load_data(model.datafile{iFile}, options.marker_chan_num);
+      [nsts, ~, ndata] = pspm_load_data(model.datafile{iSn}, options.marker_chan_num);
       if nsts == -1; warning('ID:invalid_input', 'Could not load data'); return; end
     else
-      [nsts, ~, ndata] = pspm_load_data(model.datafile{iFile}, 'marker');
+      [nsts, ~, ndata] = pspm_load_data(model.datafile{iSn}, 'marker');
       if nsts == -1; warning('ID:invalid_input', 'Could not load data'); return; end
     end
     events = ndata{1}.data;
@@ -263,7 +263,7 @@
         case 'samples'
           win = round(epochs{iSn}(iEpoch, :) * sr(datatype(k)) / sr(1));
         case 'markers'
-          win = round(events(epochs{iFile}(iEpoch, :)) * sr(datatype(k)));
+          win = round(events(epochs{iSn}(iEpoch, :)) * sr(datatype(k)));
         case 'whole'
           win = [1 numel(Y{datatype(k)})];
       end

From 4a37d1e4455aec3d15f832b3bca07ba431393d10 Mon Sep 17 00:00:00 2001
From: "teddy.chao" <dadi.zhao@ucl.ac.uk>
Date: Mon, 24 Jul 2023 14:58:18 +0100
Subject: [PATCH 43/48] Update pspm_sf.m

---
 src/pspm_sf.m | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/pspm_sf.m b/src/pspm_sf.m
index 142cb9beb..dbe0615b1 100644
--- a/src/pspm_sf.m
+++ b/src/pspm_sf.m
@@ -166,8 +166,8 @@
 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;
@@ -230,7 +230,7 @@
       [~, 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, :) = [];

From e1782e06c2547699dc945f1d6a876e7bf533b9bf Mon Sep 17 00:00:00 2001
From: "teddy.chao" <dadi.zhao@ucl.ac.uk>
Date: Mon, 24 Jul 2023 14:58:49 +0100
Subject: [PATCH 44/48] Update pspm_sf.m

---
 src/pspm_sf.m | 60 +++++++++++++++++++++++++--------------------------
 1 file changed, 30 insertions(+), 30 deletions(-)

diff --git a/src/pspm_sf.m b/src/pspm_sf.m
index 2c472ebed..763f5b1ed 100644
--- a/src/pspm_sf.m
+++ b/src/pspm_sf.m
@@ -166,8 +166,8 @@
 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;
@@ -198,15 +198,15 @@
 end
 %% 3 Get data
 missing = cell(size(model.missing));
-for iSn = 1:numel(model.datafile)
+for iFile = 1:numel(model.datafile)
   % 3.1 User output --
-  fprintf('SF analysis: %s ...', model.datafile{iSn});
+  fprintf('SF analysis: %s ...', model.datafile{iFile});
   % 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{iSn}, model.channel);
+  [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;
   model.filter.sr = sr(1);
@@ -224,51 +224,51 @@
   end
   % 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');
+  if ~isempty(model.missing{iFile})
+    [~, missing{iFile}] = pspm_get_timing('epochs', model.missing{iFile}, 'seconds');
   % 3.5.2 sort missing epochs --
-    if size(missing{iSn}, 1) > 0
-      [~, sortindx] = sort(missing{iSn}(:, 1));
-      missing{iSn} = missing{iSn}(sortindx,:);
+    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)
-        if missing{iSn}(k, 1) <= missing{iSn}(k - 1, 2)
-          missing{iSn}(k, 1) =  missing{iSn}(k - 1, 1);
-          missing{iSn}(k - 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, :) = [];
         end
       end
     end
   else
-    missing{iSn} = [];
+    missing{iFile} = [];
   end
   % 3.6 Get marker data --
   if any(strcmp(model.timeunits, {'marker', 'markers'}))
     if options.marker_chan_num
-      [nsts, ~, ndata] = pspm_load_data(model.datafile{iSn}, options.marker_chan_num);
+      [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, ~, ndata] = pspm_load_data(model.datafile{iSn}, 'marker');
+      [nsts, ~, ndata] = pspm_load_data(model.datafile{iFile}, 'marker');
       if nsts == -1; warning('ID:invalid_input', 'Could not load data'); return; end
     end
     events = ndata{1}.data;
   end
-  for iEpoch = 1:size(epochs{iSn}, 1)
+  for iEpoch = 1:size(epochs{iFile}, 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{iSn}(iEpoch, :) * sr(datatype(k)));
+          win = round(epochs{iFile}(iEpoch, :) * sr(datatype(k)));
         case 'samples'
-          win = round(epochs{iSn}(iEpoch, :) * sr(datatype(k)) / sr(1));
+          win = round(epochs{iFile}(iEpoch, :) * sr(datatype(k)) / sr(1));
         case 'markers'
-          win = round(events(epochs{iSn}(iEpoch, :)) * sr(datatype(k)));
+          win = round(events(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{iSn});
+        warning('\nEpoch %2.0f outside of file %s ...', iEpoch, model.modelfile{iFile});
       else
         % correct issues with using 'round'
         win(1) = max(win(1), 1);
@@ -276,19 +276,19 @@
       end
       % 3.6.1 collect information --
       sf.model{k}(iEpoch).modeltype = method{k};
-      sf.model{k}(iEpoch).boundaries = squeeze(epochs{iSn}(iEpoch, :));
+      sf.model{k}(iEpoch).boundaries = squeeze(epochs{iFile}(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;
-      if any(missing{iSn})
+      if any(missing{iFile})
         model.missing_data = zeros(size(escr));
-        model.missing_data((missing{iSn}(:,1)+1):(missing{iSn}(:,2)+1)) = 1;
+        model.missing_data((missing{iFile}(:,1)+1):(missing{iFile}(:,2)+1)) = 1;
       end
       % 3.6.2 do the analysis and collect results --
-      if any(missing{iSn})
+      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)));
@@ -306,16 +306,16 @@
   end
   sf.names = method(:);
   sf.infos.date = date;
-  sf.infos.file = model.modelfile{iSn};
-  sf.modelfile = model.modelfile{iSn};
+  sf.infos.file = model.modelfile{iFile};
+  sf.modelfile = model.modelfile{iFile};
   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{iSn}, 'sf');
-  outfile = model.modelfile(iSn);
+  save(model.modelfile{iFile}, 'sf');
+  outfile = model.modelfile(iFile);
   fprintf('\n');
 end
 sts = 1;

From 3cfa6490049e90b983df6e7e8cf9a7f767ca2f6a Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Mon, 24 Jul 2023 16:40:59 +0100
Subject: [PATCH 45/48] Update pspm_sf.m

---
 src/pspm_sf.m | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/pspm_sf.m b/src/pspm_sf.m
index 763f5b1ed..48c7d426f 100644
--- a/src/pspm_sf.m
+++ b/src/pspm_sf.m
@@ -285,7 +285,8 @@
       sf.model{k}(iEpoch).data = escr;
       if any(missing{iFile})
         model.missing_data = zeros(size(escr));
-        model.missing_data((missing{iFile}(:,1)+1):(missing{iFile}(:,2)+1)) = 1;
+        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})

From 3f309772e717840c735f007d85be3e8d003a2fbe Mon Sep 17 00:00:00 2001
From: "dadi.zhao" <dadi.zhao@outlook.com>
Date: Mon, 24 Jul 2023 18:32:16 +0100
Subject: [PATCH 46/48] Update pspm_split_sessions.m

---
 src/pspm_split_sessions.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pspm_split_sessions.m b/src/pspm_split_sessions.m
index 4058dc049..4c9a756aa 100644
--- a/src/pspm_split_sessions.m
+++ b/src/pspm_split_sessions.m
@@ -205,7 +205,7 @@
       newepochfile{sn} = fullfile(p_epochs, sprintf('%s_sn%02.0f%s', f_epochs, sn, ex_epochs));
     end
     % 2.4.2 Split data
-    trimoptions = struct('drop_offset_markers', 1);
+    trimoptions = struct('drop_offset_markers', 1, 'marker_chan_num', markerchannel);
     newdata = pspm_trim(struct('data', {indata}, 'infos', ininfos), ...
       options.prefix, suffix(sn), trimpoint(sn, 1:2), trimoptions);
     options.overwrite = pspm_overwrite(newdatafile{sn}, options);

From b65676222ac05cbc6055dfab096a5d427f39e16e Mon Sep 17 00:00:00 2001
From: Dominik Bach <dominik.bach@uzh.ch>
Date: Fri, 11 Aug 2023 10:25:22 +0200
Subject: [PATCH 47/48] bugfix

---
 src/pspm_cfg/pspm_cfg_sf.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pspm_cfg/pspm_cfg_sf.m b/src/pspm_cfg/pspm_cfg_sf.m
index f6b89ddf6..65d39d546 100644
--- a/src/pspm_cfg/pspm_cfg_sf.m
+++ b/src/pspm_cfg/pspm_cfg_sf.m
@@ -193,7 +193,7 @@
 mrk_chan.name    = 'Marker Channel';
 mrk_chan.tag     = 'mrk_chan';
 mrk_chan.strtype = 'i';
-mrk_chan.val     = {1};
+mrk_chan.val     = {0};
 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 ' ...

From 2cd0879de841e5340dd63aaeda2a44e523d08eb5 Mon Sep 17 00:00:00 2001
From: "teddy.chao" <dadi.zhao@ucl.ac.uk>
Date: Fri, 11 Aug 2023 17:49:48 +0100
Subject: [PATCH 48/48] Update pspm_cfg_sf.m

---
 src/pspm_cfg/pspm_cfg_sf.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pspm_cfg/pspm_cfg_sf.m b/src/pspm_cfg/pspm_cfg_sf.m
index 65d39d546..f6b89ddf6 100644
--- a/src/pspm_cfg/pspm_cfg_sf.m
+++ b/src/pspm_cfg/pspm_cfg_sf.m
@@ -193,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 ' ...