Improve pspm_trim to handle missing epoch files

src/pspm_dcm.m

@@ -363,8 +363,7 @@
 for vs = 1:numel(valid_subsessions)
   isbSn = valid_subsessions(vs);
   sbSn = subsessions(isbSn, :);
-  flanks = pspm_time2index(sbSn(2:3), sr{sbSn(1)});
-  sbSn_data = y{sbSn(1)}(flanks(1):flanks(2));
+  sbSn_data = y{sbSn(1)}((pspm_epochs2logical(sbSn(2:3), length(y{sbSn(1)}), sr{sbSn(1)})==1));
   sbs_miss = isnan(sbSn_data);
 
   if any(sbs_miss)
@@ -690,7 +689,7 @@
 
 %% 6 Invert DCM
 [sts, dcm] = pspm_dcm_inv(model, options);
-if sts < 1,
+if sts < 1
     return
 end
 

src/pspm_epochs2logical.m

@@ -15,18 +15,23 @@
 % ● Output
 %          index: [logical] index
 % ● History
-%   Introduced in PsPM 6.2
+%   Introduced in PsPM 6.1.2
 %   Written in 2024 by Dominik Bach (Uni Bonn)
 
-if nargin > 2 
-    epochs = pspm_time2index(epochs, sr, datalength);
+if nargin > 2 && sr ~= 1
+    % (1) do not let pspm_time2index account for data length, as this would
+    %     lead to wrong indices below
+    % (2) if epochs are specified in samples, then this conversion would
+    %     lead to wrong results
+    epochs = pspm_time2index(epochs, sr);
 end
 
 index = zeros(datalength, 1);
 if ~isempty(epochs)
     for k = 1:size(epochs, 1)
-        flanks = round(epochs(k,:));
-        index(flanks(1):flanks(2)) = 1;
+        flanks = epochs(k,:);
+        % ensure the epoch has duration: diff(flanks)
+        index(flanks(1):(flanks(2) - 1)) = 1;
     end
 end
 index = index(1:datalength);

src/pspm_get_timing.m

@@ -385,7 +385,7 @@
     % get epoch information from file or from input --
     if ischar(intiming)
       [sts, in] = pspm_get_timing('file', intiming);
-      if sts < 1, return; end;
+      if sts < 1, return; end
       if isfield(in, 'epochs')
         outtiming = in.epochs;
       elseif isfield(in, 'onsets')
@@ -427,8 +427,6 @@
       else
         warning('Unknown epoch definition format.');  return;
       end
-    else
-       return;
     end
 
     % check time units --

src/pspm_load_data.m

@@ -98,7 +98,7 @@
     % specify if fn is a filename
     if ~exist(fn, 'file')
       if ~isstruct(channel) % if channel is not a struct, fn must exist
-        warning('ID:nonexistent_file', 'The file fn does not exist.');
+        warning('ID:nonexistent_file', 'The file %s does not exist.', fn);
         return
       end
     else

src/pspm_options.m

@@ -451,6 +451,7 @@
     options = autofill(options, 'drop_offset_markers',    0,          '*Int'            );
     options = autofill(options, 'marker_chan_num',        'marker',   '*Int*Char'       );
     options = autofill(options, 'overwrite',              0,          1                 );
+    options = autofill(options, 'missing',               [],          '*Char'           );
   case 'write_channel'
     % 2.47 pspm_write_channel --
     options = autofill(options, 'channel',                0,          '*Int*Char'  );

src/pspm_split_sessions.m

@@ -106,27 +106,6 @@
   return;
 end
 
-% 2.2 Handle missing epochs
-if ~isempty(options.missing)
-  % makes sure the epochs are in seconds and not empty
-  [sts_get_timing, missing_time] = pspm_get_timing('epochs', options.missing, 'seconds');
-  if ~sts_get_timing
-    warning('ID:invalid_input', 'Could not load missing epochs.');
-  end
-  missingsr = 10000; % dummy sample rate, should be higher than data sampling rates (but no need to make it dynamic)
-  duration_index = round(missingsr * ininfos.duration);
-  indx = zeros(1,duration_index); % indx should be a one-dimensional array?
-  missing = pspm_time2index(missing_time, missingsr, duration_index); % convert epochs in sec to datapoints
-
-  % allow splitting empty missing epochs
-  if ~isempty(missing)
-    indx(missing(:, 1)) = 1;
-    indx(missing(:, 2)+1) = -1;
-  end
-  dp_epochs = (cumsum(indx(:)) == 1);
-  % extract fileparts for later
-  [p_epochs, f_epochs, ex_epochs] = fileparts(options.missing);
-end
 
 % 2.3 Handle markers
 
@@ -193,47 +172,27 @@
 
   % 2.4 Split files
   for sn = 1:size(trimpoint,1)
-    % 2.4.1 Determine filenames
-    [p, f, ex] = fileparts(datafile);
-    newdatafile{sn} = fullfile(p, sprintf('%s_sn%02.0f%s', f, sn, ex));
-    if ~isempty(options.missing)
-      newepochfile{sn} = fullfile(p_epochs, sprintf('%s_sn%02.0f%s', f_epochs, sn, ex_epochs));
-    end
+      % 2.4.1 Determine options & filenames
+      trimoptions = struct('drop_offset_markers', 1, 'marker_chan_num', options.marker_chan_num);
+      [p, f, ex] = fileparts(datafile);
+      newdatafile{sn} = fullfile(p, sprintf('%s_sn%02.0f%s', f, sn, ex));
+      if ~isempty(options.missing)
+        [p_epochs, f_epochs, ex_epochs] = fileparts(options.missing);
+        newepochfile{sn} = fullfile(p_epochs, sprintf('%s_sn%02.0f%s', f_epochs, sn, ex_epochs));
+        trimoptions.missing = options.missing;
+     end
     % 2.4.2 Split data
-    trimoptions = struct('drop_offset_markers', 1, 'marker_chan_num', options.marker_chan_num);
-    [tsts, newdata] = pspm_trim(struct('data', {indata}, 'infos', ininfos), ...
+    [tsts, newdata, newmissingfile] = pspm_trim(struct('data', {indata}, 'infos', ininfos), ...
       prefix{sn}, suffix{sn}, trimpoint(sn, 1:2), trimoptions);
     if tsts < 1, return; end
     options.overwrite = pspm_overwrite(newdatafile{sn}, options);
     newdata.options = options;
     pspm_load_data(newdatafile{sn}, newdata);
-    % 2.4.5 Split Epochs
+    % 2.4.3 deal with missing epoch file
     if ~isempty(options.missing)
-      dummydata{1,1}.header = struct('chantype', 'custom', ...
-        'sr', missingsr, ...
-        'units', 'unknown');
-      dummydata{1,1}.data   = dp_epochs;
-      % add marker channel so that pspm_trim has a reference
-      dummydata{2,1}      = mrkdata;
-      dummyinfos          = ininfos;
-			trimoptions_missing = trimoptions;
-			trimoptions_missing.marker_chan_num = 2;
-      [tsts, newmissing] = pspm_trim(struct('data', {dummydata}, 'infos', dummyinfos), ...
-        prefix{sn}, suffix{sn}, trimpoint(sn, 1:2), trimoptions_missing);
-      if tsts < 1, return; end
-      epochs = newmissing.data{1}.data;
-      epoch_on = 1 + strfind(epochs.', [0 1]); % Return the start points of the excluded interval
-      epoch_off = strfind(epochs.', [1 0]); % Return the end points of the excluded interval
-      if numel(epoch_off) < numel(epoch_on) % if the epochs is in the middle of 2 blocks
-        epoch_off(end+1) = numel(epochs);
-      elseif numel(epoch_on) < numel(epoch_off)
-        epoch_on = [0, epoch_on];
-      elseif (numel(epoch_on) > 0 && numel(epoch_off > 0) && epoch_off(1) < epoch_on(1))
-        epoch_on = [0, epoch_on];
-        epoch_off(end+1) = numel(epochs);
-      end
-      epochs = [epoch_on.', epoch_off.']/missingsr; % convert back to seconds
+      [sts, epochs] = pspm_get_timing('epochs', newmissingfile, 'seconds');
       save(newepochfile{sn}, 'epochs');
+      delete(newmissingfile);
     end
   end
 end

src/pspm_time2index.m

@@ -4,14 +4,14 @@
 %   to a sample index.
 % ● Format (optional arguments in []; all arguments up the last specified
 %           one need to be specified)
-%   index = pspm_time2index(time, sr [, data_length, zero_permitted, events])
+%   index = pspm_time2index(time, sr [, data_length, is_duration, events])
 % ● Arguments
 %           time: [vector or matrix] time stamps in second.
 %             sr: [numeric] sampling rate or frequency
 %    data_length: [integer] the length of data, which the index should
 %                 should not exceed
-%  zero_permitted: [0/1] whether zero is permitted (set to 0 for index and
-%                  to 1 for durations; default 0)
+%    is_duration: [0/1] whether an index or a duration is required, default
+%                 is 0
 %          events: vector of timestamps from a marker channel, will be considered if
 %                  given as input
 % ● Output
@@ -28,20 +28,23 @@
 end
 
 if nargin < 4
-    zero_permitted = 0;
+    is_duration = 0;
 else 
-    zero_permitted = 1;
+    is_duration = varargin{2};
 end
 
 if nargin > 4
     events = varargin{3};
     time = events(time);
 end
 
-index = double(int64(round(time * sr))); % round can sometimes result in non-integer values
+% 'round' can sometimes result in non-integer values
+index = double(int64(round(time * sr))); 
 
-if zero_permitted < 1
-    index(index == 0) = 1;
+% The first sample of the file corresponds to index 1 and time 0, i.e. we
+% assume the first sample was recorded in the interval [0, dt[ 
+if is_duration < 1
+    index = index + 1;
 end
 
 flag = index > ones(size(index)) * data_length;

src/pspm_trim.m

@@ -1,4 +1,4 @@
-function [sts, newdatafile] = pspm_trim(datafile, from, to, reference, options)
+function [sts, newdatafile, newepochfile] = pspm_trim(datafile, from, to, reference, options)
 % ● Description
 %   pspm_trim cuts an PsPM dataset to the limits set with the parameters 'from'
 %   and 'to' and writes it to a file with a prepended 't'
@@ -29,16 +29,23 @@
 %   │                   Default value: determined by pspm_overwrite.
 %   ├.marker_chan_num:  marker channel number.
 %   │                   if undefined or 0, first marker channel is used.
+%   ├────────.missing:  Optional name of an epoch file, e.g. containing a
+%   │                   missing epochs definition in s. This is then split
+%   │                   accordingly.
 %   └.drop_offset_markers:
-%                       if offsets are set in the reference, you might be
-%                       interested in only the data, but not in the additional
-%                       markers which are within the offset. therefore set this
+%                       if 'from' and 'to' are defined with respect to
+%                       markers, you might be interested in the data that within
+%                       extend beyond these markers but not in any additional
+%                       markers which are within this interval. Set this
 %                       option to 1 to drop markers which lie in the offset.
-%                       this is for event channels only. default is 0.
+%                       this is for event channels only. Default is 0.
 % ● Outputs
 %                  sts: status variable indicating whether function run successfully.
 %          newdatafile: a filename for the updated file (or a struct with 
 %                       fields .data and .infos if data file is a struct)
+%         newepochfile: missing epoch filename for the individual
+%                       sessions (empty if options.missing not specified)
+
 % ● Version
 %   Introduced in PsPM 3.0
 %   Written in 2008-2015 by Dominik R Bach (Wellcome Trust Centre for Neuroimaging)
@@ -52,6 +59,7 @@
 end
 sts = -1;
 newdatafile = [];
+newepochfile = [];
 
 % 1.2 Verify the number of input arguments
 switch nargin
@@ -263,65 +271,73 @@
         sta_offset = 0;
     end
 end
+sta_time = sta_p + sta_offset;
 if (sto_p + sto_offset) > infos.duration
     warning('ID:marker_out_of_range', ['\nEnd point (%.2f s) outside ', ...
         'file, no trimming at end.'], (sto_p + sto_offset));
-    if (sto_p > infos.duration)
-        sto_p = infos.duration;
-        sto_offset = 0;
-    else
-        sto_offset = infos.duration - sto_p;
-    end
+    % adjustment of the end point is being taken care of in pspm_epochs2logical
+    sto_time = infos.duration;
+else
+    sto_time = sto_p + sto_offset;
 end
+
 % 2.5 Trim file
 for k = 1:numel(data)
     if ~strcmpi(data{k}.header.units, 'events') % waveform channels
-        % set start point (`ceil` for protect against having duration < data*sr,
-        % the `+1` is here because of matlabs convention to start indices from 1)
-        newstartpoint = ceil((sta_p + sta_offset) * data{k}.header.sr)+1;
-        if newstartpoint == 0
-            newstartpoint = 1;
-        end
-        % set end point
-        newendpoint = floor((sto_p + sto_offset) * data{k}.header.sr);
-        if newendpoint > numel(data{k}.data), ...
-                newendpoint = numel(data{k}.data);
-        end
-        % trim data
-        data{k}.data=data{k}.data(newstartpoint:newendpoint);
+        index = pspm_epochs2logical([sta_time, sto_p + sto_offset], ...
+                                     numel(data{k}.data), ...
+                                     data{k}.header.sr);
+        data{k}.data=data{k}.data(find(index));
     else % event channels
         if options.drop_offset_markers
-            newendpoint = sto_p;
             newstartpoint = sta_p;
+            newendpoint   = sto_p;
         else
-            newendpoint = sto_p + sto_offset;
-            newstartpoint = sta_p + sta_offset;
+            newstartpoint = sta_time;
+            newendpoint   = sto_time;
         end
+        newendpoint = min([newendpoint, sto_time]);
+        remove_early = data{k}.data < newstartpoint;
         remove_late = data{k}.data > newendpoint;
-        data{k}.data(remove_late) = [];
-        data{k}.data = data{k}.data - newstartpoint;
-        remove_early = data{k}.data < 0;
-        data{k}.data(remove_early) = [];
-
-        % move to match data if offset markers should be dropped
-        if options.drop_offset_markers
-            data{k}.data = data{k}.data - sta_offset;
-        end
+        remove_index = any([remove_early, remove_late], 2);
+        data{k}.data(remove_index) = [];
+        data{k}.data = data{k}.data - sta_time;
+
         if isfield(data{k}, 'markerinfo')
             % also trim marker info if available
-            data{k}.markerinfo.value(remove_late) = [];
-            data{k}.markerinfo.name(remove_late) = [];
-
-            data{k}.markerinfo.value(remove_early) = [];
-            data{k}.markerinfo.name(remove_early) = [];
+            data{k}.markerinfo.value(remove_index) = [];
+            data{k}.markerinfo.name(remove_index) = [];
         end
     end
     % save new file
-    infos.duration = (sto_p + sto_offset) - (sta_p + sta_offset);
+    infos.duration = sto_time - sta_time;
     infos.trimdate = date;
-    infos.trimpoints = [(sta_p + sta_offset) (sto_p + sto_offset)];
+    infos.trimpoints = [sta_time, sto_time];
 end
 clear savedata
+
+% handle optional missing data file
+if ~isempty(options.missing)
+    [lsts, epochs] = pspm_get_timing('epochs', options.missing, 'seconds');
+    if lsts < 1, return; end
+    index = epochs(:, 2) < sta_time | ...
+            epochs(:, 1) > sto_time | ...
+            epochs(:, 1) > infos.duration;
+    epochs(index, :) = [];
+    epochs = epochs - sta_time;
+    if ~isempty(epochs)
+        epochs(1, 1) = max([0, epochs(1, 1)]);
+        epochs(end, 2) = min([infos.duration, epochs(end, 2)]);
+    end
+    lsts = pspm_get_timing('epochs', epochs, 'seconds');
+    if lsts < 1, return; end
+    [pth, fn, ext] = fileparts(options.missing);
+    newepochfile = fullfile(pth, ['t', fn, ext]);
+    save(newepochfile, 'epochs');
+end
+
+
+
 % 2.6 Save data
 savedata.data = data;
 savedata.infos = infos;