function [data] = ft_datatype_raw(data, varargin)

% FT_DATATYPE_RAW describes the FieldTrip MATLAB structure for raw data
%
% The raw datatype represents sensor-level time-domain data typically
% obtained after calling FT_DEFINETRIAL and FT_PREPROCESSING. It contains
% one or multiple segments of data, each represented as Nchan X Ntime
% arrays.
%
% An example of a raw data structure with 151 MEG channels is
%
%          label: {151x1 cell}      the channel labels represented as a cell-array of strings
%           time: {1x266 cell}      the time axis [1*Ntime double] per trial
%          trial: {1x266 cell}      the numeric data as a cell array, with a matrix of [151*Ntime double] per trial
%     sampleinfo: [266x2 double]    the begin and endsample of each trial relative to the recording on disk
%      trialinfo: [266x1 double]    optional trigger or condition codes for each trial
%            hdr: [1x1 struct]      the full header information of the original dataset on disk
%           grad: [1x1 struct]      information about the sensor array (for EEG it is called elec)
%            cfg: [1x1 struct]      the configuration used by the function that generated this data structure
%
% Required fields:
%   - time, trial, label
%
% Optional fields:
%   - sampleinfo, trialinfo, grad, elec, opto, hdr, cfg
%
% Deprecated fields:
%   - fsample
%
% Obsoleted fields:
%   - offset
%
% Revision history:
%
% (2011/latest) The description of the sensors has changed, see FT_DATATYPE_SENS
% for further information.
%
% (2010v2) The trialdef field has been replaced by the sampleinfo and
% trialinfo fields. The sampleinfo corresponds to trl(:,1:2), the trialinfo
% to trl(4:end).
%
% (2010v1) In 2010/Q3 it shortly contained the trialdef field which was a copy
% of the trial definition (trl) is generated by FT_DEFINETRIAL.
%
% (2007) It used to contain the offset field, which corresponds to trl(:,3).
% Since the offset field is redundant with the time axis, the offset field is
% from now on not present any more. It can be recreated if needed.
%
% (2003) The initial version was defined
%
% See also FT_DATATYPE, FT_DATATYPE_COMP, FT_DATATYPE_TIMELOCK, FT_DATATYPE_FREQ,
% FT_DATATYPE_SPIKE, FT_DATATYPE_SENS

% Copyright (C) 2011, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% get the optional input arguments, which should be specified as key-value pairs
version       = ft_getopt(varargin, 'version', 'latest');
hassampleinfo = ft_getopt(varargin, 'hassampleinfo', 'ifmakessense'); % can be yes/no/ifmakessense
hastrialinfo  = ft_getopt(varargin, 'hastrialinfo',  'ifmakessense'); % can be yes/no/ifmakessense
allowemptytrials = ft_getopt(varargin, 'allowemptytrials', 'yes');

% do some sanity checks
assert(isfield(data, 'trial') && isfield(data, 'time') && isfield(data, 'label'), 'inconsistent raw data structure, some field is missing');
assert(length(data.trial)==length(data.time), 'inconsistent number of trials in raw data structure');
for i=1:length(data.trial)
  assert(size(data.trial{i},2)==length(data.time{i}), 'inconsistent number of samples in trial %d', i);
  assert(size(data.trial{i},1)==length(data.label), 'inconsistent number of channels in trial %d', i);
end
assert(length(unique(data.label))==length(data.label), 'channel labels must be unique');

% convert it into true/false
if isequal(hassampleinfo, 'ifmakessense')
  hassampleinfo = makessense(data, 'sampleinfo');
else
  hassampleinfo = istrue(hassampleinfo);
end
if isequal(hastrialinfo, 'ifmakessense')
  hastrialinfo = makessense(data, 'trialinfo');
else
  hastrialinfo = istrue(hastrialinfo);
end
allowemptytrials = istrue(allowemptytrials);

if strcmp(version, 'latest')
  version = '2024';
end

if isempty(data)
  return;
end

switch version
  case '2024'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if ~allowemptytrials
      remove = false(numel(data.time),1);
      for i=1:numel(data.time)
        if isempty(data.time{i})
          remove(i) = true;
        end
      end
      data.time  = data.time(~remove);
      data.trial = data.trial(~remove);
      if isfield(data, 'sampleinfo')
        data.sampleinfo = data.sampleinfo(~remove,:);
      end
      if isfield(data, 'trialinfo')
        data.trialinfo = data.trialinfo(~remove,:);
      end
    end

    % ensure that the sensor structures are up to date
    if isfield(data, 'grad')
      data.grad = ft_datatype_sens(data.grad);
    end
    if isfield(data, 'elec')
      data.elec = ft_datatype_sens(data.elec);
    end
    if isfield(data, 'opto')
      data.opto = ft_datatype_sens(data.opto);
    end
    
    if ~isfield(data, 'fsample')
      for i=1:length(data.time)
        if length(data.time{i})>1
          data.fsample = 1/mean(diff(data.time{i}));
          break
        else
          data.fsample = nan;
        end
      end
      if isnan(data.fsample)
        ft_warning('cannot determine sampling frequency');
      end
    end
    
    if isfield(data, 'offset')
      data = rmfield(data, 'offset');
    end
    
    % the trialdef field should be renamed into sampleinfo
    if isfield(data, 'trialdef')
      data.sampleinfo = data.trialdef;
      data = rmfield(data, 'trialdef');
    end
    
    if (hassampleinfo && ~isfield(data, 'sampleinfo')) || (hastrialinfo && ~isfield(data, 'trialinfo'))
      % try to reconstruct the sampleinfo and trialinfo
      data = fixsampleinfo(data);
    end
    
    if ~hassampleinfo && isfield(data, 'sampleinfo')
      data = rmfield(data, 'sampleinfo');
    end
    
    if ~hastrialinfo && isfield(data, 'trialinfo')
      data = rmfield(data, 'trialinfo');
    end
    
    if isfield(data, 'sampleinfo') && istable(data.sampleinfo)
      % the sampleinfo contains two columns with the begsample and endsample and can always be represented as a numeric array
      % the trialinfo can be either a numeric array or a table
      data.sampleinfo = table2array(data.sampleinfo);
    end
    
  case '2011'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % ensure that the sensor structures are up to date
    if isfield(data, 'grad')
      data.grad = ft_datatype_sens(data.grad);
    end
    if isfield(data, 'elec')
      data.elec = ft_datatype_sens(data.elec);
    end
    if isfield(data, 'opto')
      data.opto = ft_datatype_sens(data.opto);
    end
    
    if ~isfield(data, 'fsample')
      for i=1:length(data.time)
        if length(data.time{i})>1
          data.fsample = 1/mean(diff(data.time{i}));
          break
        else
          data.fsample = nan;
        end
      end
      if isnan(data.fsample)
        ft_warning('cannot determine sampling frequency');
      end
    end
    
    if isfield(data, 'offset')
      data = rmfield(data, 'offset');
    end
    
    % the trialdef field should be renamed into sampleinfo
    if isfield(data, 'trialdef')
      data.sampleinfo = data.trialdef;
      data = rmfield(data, 'trialdef');
    end
    
    if (hassampleinfo && ~isfield(data, 'sampleinfo')) || (hastrialinfo && ~isfield(data, 'trialinfo'))
      % try to reconstruct the sampleinfo and trialinfo
      data = fixsampleinfo(data);
    end
    
    if ~hassampleinfo && isfield(data, 'sampleinfo')
      data = rmfield(data, 'sampleinfo');
    end
    
    if ~hastrialinfo && isfield(data, 'trialinfo')
      data = rmfield(data, 'trialinfo');
    end
    
    if isfield(data, 'sampleinfo') && istable(data.sampleinfo)
      % the sampleinfo contains two columns with the begsample and endsample and can always be represented as a numeric array
      % the trialinfo can be either a numeric array or a table
      data.sampleinfo = table2array(data.sampleinfo);
    end
    
  case '2010v2'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if ~isfield(data, 'fsample')
      data.fsample = 1/mean(diff(data.time{1}));
    end
    
    if isfield(data, 'offset')
      data = rmfield(data, 'offset');
    end
    
    % the trialdef field should be renamed into sampleinfo
    if isfield(data, 'trialdef')
      data.sampleinfo = data.trialdef;
      data = rmfield(data, 'trialdef');
    end
    
    if (hassampleinfo && ~isfield(data, 'sampleinfo')) || (hastrialinfo && ~isfield(data, 'trialinfo'))
      % try to reconstruct the sampleinfo and trialinfo
      data = fixsampleinfo(data);
    end
    
    if ~hassampleinfo && isfield(data, 'sampleinfo')
      data = rmfield(data, 'sampleinfo');
    end
    
    if ~hastrialinfo && isfield(data, 'trialinfo')
      data = rmfield(data, 'trialinfo');
    end
    
  case {'2010v1' '2010'}
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if ~isfield(data, 'fsample')
      data.fsample = 1/mean(diff(data.time{1}));
    end
    
    if isfield(data, 'offset')
      data = rmfield(data, 'offset');
    end
    
    if ~isfield(data, 'trialdef') && hascfg
      % try to find it in the nested configuration history
      data.trialdef = ft_findcfg(data.cfg, 'trl');
    end
    
  case '2007'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if ~isfield(data, 'fsample')
      data.fsample = 1/mean(diff(data.time{1}));
    end
    
    if isfield(data, 'offset')
      data = rmfield(data, 'offset');
    end
    
  case '2003'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if ~isfield(data, 'fsample')
      data.fsample = 1/mean(diff(data.time{1}));
    end
    
    if ~isfield(data, 'offset')
      data.offset = zeros(length(data.time),1);
      for i=1:length(data.time);
        data.offset(i) = round(data.time{i}(1)*data.fsample);
      end
    end
    
  otherwise
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ft_error('unsupported version "%s" for raw datatype', version);
end


% Numerical inaccuracies in the binary representations of floating point
% values may accumulate. The following code corrects for small inaccuracies
% in the time axes of the trials. See http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=1390
data = fixtimeaxes(data);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function data = fixtimeaxes(data)

if ~isfield(data, 'fsample')
  fsample = 1/mean(diff(data.time{1}));
else
  fsample = data.fsample;
end

begtime   = nan(1, length(data.time));
endtime   = nan(1, length(data.time));
numsample = zeros(1, length(data.time));
for i=1:length(data.time)
  if ~isempty(data.time{i})
    begtime(i) = data.time{i}(1);
    endtime(i) = data.time{i}(end);
  end
  numsample(i) = length(data.time{i});
end

% compute the differences over trials and the tolerance
tolerance     = 0.01*(1/fsample);
begdifference = abs(begtime-begtime(1));
enddifference = abs(endtime-endtime(1));

% check whether begin and/or end are identical, or close to identical
begidentical  = all(begdifference==0);
endidentical  = all(enddifference==0);
begsimilar    = all(begdifference < tolerance);
endsimilar    = all(enddifference < tolerance);

% Compute the offset of each trial relative to the first trial, and express
% that in samples. Non-integer numbers indicate that there is a slight skew
% in the time over trials. This works in case of variable length trials.
offset = fsample * (begtime-begtime(1));
skew   = abs(offset - round(offset));

% try to determine all cases where a correction is needed
% note that this does not yet address all possible cases where a fix might be needed
needfix = false;
needfix = needfix || ~begidentical && begsimilar;
needfix = needfix || ~endidentical && endsimilar;
needfix = needfix || ~all(skew==0) && all(skew<0.01);

% if the skew is less than 1% it will be corrected
if needfix
  ft_warning('correcting numerical inaccuracy in the time axes');
  for i=1:length(data.time)
    % reconstruct the time axis of each trial, using the begin latency of
    % the first trial and the integer offset in samples of each trial
    data.time{i} = begtime(1) + ((1:numsample(i)) - 1 + round(offset(i)))/fsample;
  end
end
