% JOIN multiple AOs into a single AO.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: JOIN multiple AOs into a single AO.
%              If any two AOs overlap, then the values from the first appear
%              in the output AO. 
% 
%              Note: If the input AOs are of type 'tsdata', then they will
%              be sorted in ascending order according the t0 of each
%              object.
% 
%
% CALL:        bs = join(a1,a2,a3,...,pl)
%              bs = join(as,pl)
%              bs = as.join(pl)
%
% INPUTS:      aN   - input analysis objects
%              as   - input analysis objects array
%              pl   - input parameter list
%
% OUTPUTS:     b    - output analysis object
%
% <a href="matlab:utils.helper.displayMethodInfo('ao', 'join')">Parameters Description</a>
%
% REMARK:      Input AOs should be of the same type; if not, only AOs of the
%              type of the first input AO will be joined together to produce
%              the output.
%
% VERSION:     $Id: join.m,v 1.60 2011/09/02 11:11:51 ingo Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% PARAMETERS:  'zerofill' - Fills with zeros the gaps between the data
%                           points of the subsequent aos. [Default: 'no']
%             'sameT0'    - Does not recalculate t0 but uses the common
%                           one. [Default: 'no']
%                           Note: the t0 among different objects must be the same!
%

function varargout = join(varargin)
  
  % Check if this is a call for parameters
  if utils.helper.isinfocall(varargin{:})
    varargout{1} = getInfo(varargin{3});
    return
  end
  
  import utils.const.*
  utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
  
  % Collect input variable names
  in_names = cell(size(varargin));
  for ii = 1:nargin,in_names{ii} = inputname(ii);end
  
  % Collect all AOs
  [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
  pl             = utils.helper.collect_objects(varargin(:), 'plist', in_names);
  
  % Combine plists
  pl = parse(pl, getDefaultPlist);
  

  
  %----------------------------------------------
  % Get data type from the first AO
  dtype = class(as(1).data);
  
  % Sort the input AOs by t0, if applicable
  if strcmp(dtype, 'tsdata')
    times = as.t0.double;
    [~, idx] = sort(times);
    as = as(idx);
  end
  
  %----------------------------------------------
  % Go through each AO and collect the data of type 'dtype'
  
  histin = [];
  xo     = [];
  yo     = [];
  dxo    = [];
  dyo    = [];
  enbw0  = [];
  fs     = -1;
  aname  = '';
  adescr = '';
  plotinfo = [];
  if as(1).data.isprop('xunits')
    xunitsSimple = simplify(as(1).data.xunits);
    xunits       = as(1).data.xunits;
  end
  yunitsSimple = simplify(as(1).data.yunits);
  yunits       = as(1).data.yunits;
  
  % Get the tolerance for considering fs equal
  fstol = find(pl, 'fstol');

  % Compute time offset for tsdata objects to avoid rounding errors later
  minT0milli = getMinT0(as);
  
  % loop over AOs
  for jj=1:numel(as)
    % Only get the data type we want
    if isa(as(jj).data, dtype)
      switch lower(dtype)
        case 'tsdata'
          % here we concatonate time-series
          t0 = (as(jj).data.t0.utc_epoch_milli - minT0milli)/1000;
          % make proper time vector
          x = as(jj).x + t0;
          
          % only add samples past the end of existing (first loop)
          if isempty(xo)
            yo = as(jj).y;
            xo = x;
            if numel(as(jj).dx) == 0
              dxo = zeros(numel(as(jj).x),1);
            elseif numel(as(jj).dx) == 1
              dxo = ones(numel(as(jj).x),1) .* as(jj).dx;
            else
              dxo = as(jj).dx;
            end
            if numel(as(jj).dy) == 0
              dyo = zeros(numel(as(jj).y),1);
            elseif numel(as(jj).dy) == 1
              dyo = ones(numel(as(jj).y),1) .* as(jj).dy;
            else
              dyo = as(jj).dy;
            end
          else
            idxPost = find(x > max(xo));
            idxPre  = find(x < min(xo));
            
            %%%%%%%%%%   Fill the gaps with zeros   %%%%%%%%%%
            zerofill = utils.prog.yes2true(find(pl, 'zerofill'));
            
            if zerofill
              % Check if there is a gap between the x-values and the pre-values.
              if ~isempty(idxPre)
                interStart = x(idxPre(end));
                interEnd   = xo(1);
                nsecsPre2no = interEnd - interStart;
                
                % The gap must be larger than 1/fs in order to
                % fill the gap with zeros
                if nsecsPre2no > 1/fs
                  x_interPre = linspace(interStart+1/fs, interEnd-1/fs, nsecsPre2no*fs-2*1/fs).';
                  y_interPre = zeros(length(x_interPre), 1);
                else
                  x_interPre = [];
                  y_interPre = [];
                end
              else
                x_interPre = [];
                y_interPre = [];
              end
              
              % Check if there is a gap between the x-values and the post-values.
              if ~isempty(idxPost)
                interStart   = xo(end);
                interEnd     = x(idxPost(1));
                nsecsPost2no = interEnd - interStart;
                
                % The gap must be larger than 1/fs in order to
                % fill the gap with zeros
                if nsecsPost2no > 1/fs
                  x_interPost = linspace(interStart+1/fs, interEnd-1/fs, nsecsPost2no*fs-1/fs).';
                  y_interPost = zeros(length(x_interPost), 1);
                else
                  x_interPost = [];
                  y_interPost = [];
                end
              else
                x_interPost = [];
                y_interPost = [];
              end
              
            else
              %%%%%%%%%%   Don't fill the gaps with zeros   %%%%%%%%%%
              x_interPre  = [];
              y_interPre  = [];
              x_interPost = [];
              y_interPost = [];
            end
            xo = [x(idxPre); x_interPre; xo; x_interPost; x(idxPost)];
            yo = [as(jj).data.getY(idxPre); y_interPre; yo; y_interPost; as(jj).data.getY(idxPost)];
            
            %%% Collect errors
            if numel(as(jj).dx) == 0
              dx = zeros(numel(as(jj).x),1);
            elseif numel(as(jj).dx) == 1
              dx = ones(numel(as(jj).x),1) .* as(jj).dx;
            else
              dx = as(jj).dx;
            end
            
            if numel(as(jj).dy) == 0
              dy = zeros(numel(as(jj).y),1);
            elseif numel(as(jj).dy) == 1
              dy = ones(numel(as(jj).y),1) .* as(jj).dy;
            else
              dy = [dyo; as(jj).dy];
            end
            
            x_interPre = zeros(numel(x_interPre),1);
            y_interPre = zeros(numel(y_interPre),1);
            x_interPost = zeros(numel(x_interPost),1);
            y_interPost = zeros(numel(y_interPost),1);
            dxo = [dx(idxPre); x_interPre; dxo; x_interPost; dx(idxPost)];
            dyo = [dy(idxPre); y_interPre; dyo; y_interPost; dy(idxPost)];
            
          end
          
          % check fs
          if (fs > 0) && (abs(as(jj).fs - fs) > fstol*fs)
            error('### Data has different sample rates');
          end
          % store fs
          fs = as(jj).fs;
          
          % check xunits
          if ~eq(xunitsSimple, simplify(as(jj).xunits))
            error('### The x-units of the analysis objects are not the same %s <-> %s', char(xunits), char(as(jj).xunits));
          end
          % check yunits
          if ~eq(yunitsSimple, simplify(as(jj).yunits))
            error('### The y-units of the analysis objects are not the same %s <-> %s', char(yunits), char(as(jj).yunits));
          end
          
          % store T0
          T0s(jj) = as(jj).t0;
          
        case 'fsdata'
          
          %%% Collect all fsdata samples
          if isempty(xo)
            idxBefore = 1:numel(as(jj).x);
            idxAfter  = [];
          else
            idxBefore = find(as(jj).x < xo(1));
            idxAfter  = find(as(jj).x > xo(end));
          end
          xo = [as(jj).x(idxBefore); xo; as(jj).x(idxAfter)];
          yo = [as(jj).y(idxBefore); yo; as(jj).y(idxAfter)];

          %%% Collect all errors
          % dx
          if numel(as(jj).dx) == 0
            dx = zeros(numel(as(jj).x),1);
          elseif numel(as(jj).dx) == 1
            dx = ones(numel(as(jj).x),1) .* as(jj).dx;
          else
            dx = as(jj).dx;
          end
          dxo = [dx(idxBefore); dxo; dx(idxAfter)];
          % dy
          if numel(as(jj).dy) == 0
            dy = zeros(numel(as(jj).y),1);
          elseif numel(as(jj).dy) == 1
            dy = ones(numel(as(jj).y),1) .* as(jj).dy;
          else
            dy = as(jj).dy;
          end
          dyo = [dy(idxBefore); dyo; dy(idxAfter)];
          % enbw
          if numel(as(jj).data.enbw) == 0
            enbw = NaN(numel(as(jj).y),1);
          elseif numel(as(jj).data.enbw) == 1
            enbw = ones(numel(as(jj).y),1) .* as(jj).data.enbw;
          else
            enbw = as(jj).data.enbw;
          end
          enbw0 = [enbw(idxBefore); enbw0; enbw(idxAfter)];

          % check fs
          if (fs > 0) && (abs(as(jj).fs - fs) > fstol*fs)
            error('### Data has different sample rates');
          end
          % store fs
          fs = as(jj).fs;
          
          % check xunits
          if ~eq(xunitsSimple, simplify(as(jj).xunits))
            error('### The x-units of the analysis objects are not the same %s <-> %s', char(xunits), char(as(jj).xunits));
          end
          % check yunits
          if ~eq(yunitsSimple, simplify(as(jj).yunits))
            error('### The y-units of the analysis objects are not the same %s <-> %s', char(yunits), char(as(jj).yunits));
          end
          
        case 'xydata'
          xo = [xo; as(jj).x];
          yo = [yo; as(jj).y];
          if numel(as(jj).dx) == 0
            dx = zeros(numel(as(jj).x),1);
          elseif numel(as(jj).dx) == 1
            dx = ones(numel(as(jj).x),1) .* as(jj).dx;
          else
            dx = as(jj).dx;
          end
          dxo = [dxo; dx];
          if numel(as(jj).dy) == 0
            dy = zeros(numel(as(jj).y),1);
          elseif numel(as(jj).dy) == 1
            dy = ones(numel(as(jj).y),1) .* as(jj).dy;
          else
            dy = as(jj).dy;
          end
          dyo = [dyo; dy];
          
          % check xunits
          if ~eq(xunitsSimple, simplify(as(jj).xunits))
            error('### The x-units of the analysis objects are not the same %s <-> %s', char(xunits), char(as(jj).xunits));
          end
          % check yunits
          if ~eq(yunitsSimple, simplify(as(jj).yunits))
            error('### The y-units of the analysis objects are not the same %s <-> %s', char(yunits), char(as(jj).yunits));
          end
          
        case 'cdata'
          try
            yo = [yo; as(jj).y];
            if numel(as(jj).dy) == 0
              dy = zeros(numel(as(jj).y),1);
            elseif numel(as(jj).dy) == 1
              dy = ones(numel(as(jj).y),1) .* as(jj).dy;
            else
              dy = as(jj).dy;
            end
            dyo = [dyo; dy];
          catch E
            disp(E.message)
            error('### It is not possible to join the data or error because they have different dimensions.');
          end
          
          % check yunits
          if ~eq(yunitsSimple, simplify(as(jj).yunits))
            error('### The y-units of the analysis objects are not the same %s <-> %s', char(yunits), char(as(jj).yunits));
          end
          
        otherwise
          error('### Unknown data type');
      end
      % Collect this input history
      histin = [histin as(jj).hist];
      % Collect the 'plotinfo'
      if ~isempty(as(jj).plotinfo)
        plotinfo = combine(plotinfo, as(jj).plotinfo);
      end
      % Collect the descriptions
      adescr = strcat(adescr, as(jj).description);
      % Collect names, invars
      if ~isempty(aname)
        if ~strcmp(aname, as(jj).name)
          aname = [aname ',' as(jj).name];
        end
      else
        aname = as(jj).name;
      end
    else
      warning('!!! Ignoring AO input with data type %s', dtype);
    end
  end
  
  %----------------------------------------------
  % Now sort output vectors
  if ~isempty(xo)
    [xos, idx] = sort(xo);
    yos = yo(idx);
    dxos = dxo(idx);
    dyos = dyo(idx);
  else
    xos = xo;
    yos = yo;
    dxos = dxo;
    dyos = dyo;
  end
  
  if all(dxos == 0)
    dxos = [];
  elseif all(diff(dxos) == 0)
    dxos = dxos(1);
  end
  if all(dyos == 0)
    dyos = [];
  elseif all(diff(dyos) == 0)
    dyos = dyos(1);
  end
  
  % Keep the data shape if the input AO
  if size(as(1).data.y,1) == 1
    xos = xos.';
    yos = yos.';
  end
  
  %%% Build output data object
  switch lower(dtype)
    case 'tsdata'
      sameT0 = utils.prog.yes2true(find(pl, 'sameT0'));
      
%       % Check that all objects share the t0
%       for kk = 1:numel(T0s)
%         if ne(T0s(kk), T0s(1))
%           sameT0 = false;
%         end
%       end
%       if sameT0
%         data = tsdata(xos, yos, fs);
%         data.setT0(minT0milli/1000);
%       else
        % get t0
%         toffset  = xos(1);
%         xos = xos - toffset;
        data = tsdata(xos, yos, fs);
        data.setT0((minT0milli/1000));
%         data.setToffset(toffset*1000);
%       end
      data.setDx(dxos);
      data.setDy(dyos);
      data.setXunits(xunits);
      data.setYunits(yunits);
      data.collapseX;
    case 'fsdata'
      if all(isnan(enbw0))
        enbw0 = [];
      elseif all(diff(enbw0) == 0)
        enbw0 = enbw0(1);
      end
      data = fsdata(xos, yos);
      data.setDx(dxos);
      data.setDy(dyos);
      data.setEnbw(enbw0);
      data.setFs(fs);
      data.setXunits(xunits);
      data.setYunits(yunits);
    case 'xydata'
      data = xydata(xos, yos);
      data.setDx(dxos);
      data.setDy(dyos);
      data.setXunits(xunits);
      data.setYunits(yunits);
    case 'cdata'
      data = cdata(yos);
      data.setDy(dyos);
      data.setYunits(yunits);
  end
  
  %----------------------------------------------
  % Build output AO
  if nargout == 0
    a = as(1);
    a.data = data;
  else
    a = ao(data);
  end
  % Set name
  a.name = aname;
  % Set description
  a.description = adescr;
  % Set plotinfo
  a.plotinfo = plotinfo;
  % Add history
  a.addHistory(getInfo('None'), pl, ao_invars, histin);
  
  %----------------------------------------------
  % set output
  varargout{1} = a;
end

%--------------------------------------------------------------------------
% Get Info Object
%--------------------------------------------------------------------------
function ii = getInfo(varargin)
  if nargin == 1 && strcmpi(varargin{1}, 'None')
    sets = {};
    pl   = [];
  else
    sets = {'Default'};
    pl   = getDefaultPlist;
  end
  % Build info object
  ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.helper, '$Id: join.m,v 1.60 2011/09/02 11:11:51 ingo Exp $', sets, pl);
  ii.setArgsmin(2);
end

%--------------------------------------------------------------------------
% Get Default Plist
%--------------------------------------------------------------------------
function plout = getDefaultPlist()
  persistent pl;  
  if exist('pl', 'var')==0 || isempty(pl)
    pl = buildplist();
  end
  plout = pl;  
end

function pl = buildplist()
  
  pl = plist();
  
  % Zero fill
  p = param({'zerofill','Fills with zeros the gaps between the data points of the subsequent aos.'}, paramValue.YES_NO);
  p.val.setValIndex(2);
  pl.append(p);
  
  % Same T0
  p = param({'sameT0', ['Does not recalculate t0 but uses the common one.<br>', ...
    'Note: the t0 among different objects must be the same!']}, paramValue.YES_NO);
  p.val.setValIndex(2);
  pl.append(p);
  
  % fstol
  p = param({'fstol', ['Relative tolerance between sampling frequency of different objects.<br>', ...
    'Jitter in the sampling frequency by less than this amount will be neglected.<br>', ...
    'If the difference is more than the set value, an error will occur.']}, paramValue.DOUBLE_VALUE(1e-6));
  pl.append(p);
end


%--------------------------------------------------------------------------
% Get Offset of this set of time-vectors
%--------------------------------------------------------------------------
function Toff = getMinT0(as)
  Toff = min(double(as.t0))*1e3;
end