function tscat = concat_timeseries(varargin)

% Concatenate multiple timeseries structures (e.g., from different runs)
% into one structure. All structures must have the same number of ROIs, and
% ROIs must be the same order. Each structure does not have to contain the
% same set of participants, nor must they be in the same order. Several
% options are available for normalizing each timeseries before
% concatenation, including: no normalization, mean-centering (default),
% z-scoring. We recommend against no normalization, as this may have a
% strong (and differential by participant) impact on association measures.
% Arguments can be input in any order. 
% 
% Example use:
%             tscat = concat_timeseries(ts1,2,ts2,ts3);
%
%             Where ts1, ts2, and ts3 are structures and 2 specifies
%             z-scoring each timeseries before concatenation. 
%
%
% Inputs: ts,    structures with timeseries data. Each structure must have
%                timeseries data as a field labeled ts, along with cell
%                array of participant IDs as a field labeled subs
%         flag,  specify what normalization should occur before
%                concatenation: 
%                              0 = do nothing
%                              1 = mean center each timeseries individually
%                                  before concatenation (to avoid mean
%                                  signal differences between runs driving
%                                  association measures) (default)
%                              2 = zscore each timeseries individually
%                                  before concatenation (to additionally
%                                  avoid higher variance timeseries having
%                                  a greater impact on association
%                                  measures)
%
% Output: tscat, structure with concatenated timeseries. Be sure to check
%                the order of participants, as this can differ from the
%                inputs.
%
% Author: Jeffrey M. Spielberg (jspielb2@gmail.com)
% Version: 01.17.17
%
% WARNING: This is a beta version. There are no known bugs, but only
% limited  testing has been perfomed. This software comes with no warranty
% (even the implied warranty of merchantability or fitness for a particular
% purpose). Therefore, USE AT YOUR OWN RISK!!!
%
% Copyleft 2014-2017. Software can be modified and redistributed, but 
% modifed, redistributed versions must have the same rights

flag  = 1;                         % Default to mean-centering before concatenation
tss   = cell(2,1);                 % Preallocate space
count = 1;                         % Start counter
for ca = 1:nargin                  % Parse input arguments
    if isstruct(varargin{ca})      % If current argument is a structure
        tss{count} = varargin{ca}; % Set as timeseries
        count      = count+1;      % Increment timeseries
    elseif isnumeric(varargin{ca}) % If current argument is an integer
        flag = varargin{ca};       % Set as option for normalization
    end
end

% Check whether each input has the same number of ROIs
nROIs = size(tss{1}.ts{1},1);
for cts = 2:size(tss,1)
    if nROIs~=size(tss{cts}.ts{1},1)
        error('Timeseries must have the same number of ROIs')
    end
end

% Find ROI labels if present
for cts = 1:size(tss,1)
    if isfield(tss{cts},'ROI_labels')
        tscat.ROI_labels = tss{cts}.ROI_labels;
    end
end

for csub = 1:length(tss{1}.subs)                                                                                                  % For each participant in first input
    tscat.subs{csub,1} = tss{1}.subs{csub};                                                                                       % Set current identifier
    switch flag
        case 0
            tscat.ts{csub,1} = tss{1}.ts{csub};                                                                                   % Perform no normalization
        case 1 
            for cROI = 1:nROIs                                                                                                    % For each ROI
                tscat.ts{csub,1}(cROI,:) = tss{1}.ts{csub}(cROI,:)-mean(tss{1}.ts{csub}(cROI,:));                                 % Mean center each timeseries before concatenation
            end
        case 2
            for cROI = 1:nROIs                                                                                                    % For each ROI
                tscat.ts{csub,1}(cROI,:) = (tss{1}.ts{csub}(cROI,:)-mean(tss{1}.ts{csub}(cROI,:)))./std(tss{1}.ts{csub}(cROI,:)); % z-score each timeseries before concatenation
            end
    end
end

for cts = 2:size(tss,1)                                                                                                                                          % For each timeseries input after the first
    for csub = 1:length(tss{cts}.subs)                                                                                                                           % For each participant in current structure
        pm = ismember(tscat.subs,tss{cts}.subs{csub});                                                                                                           % Determine whether current participant is present in currently concatenated data
        if any(pm)
            switch flag
                case 0
                    tscat.ts{pm,1} = [tscat.ts{pm,1},tss{cts}.ts{csub}];                                                                                         % Perform no normalization
                case 1
                    ntp1 = size(tscat.ts{pm,1},2);                                                                                                               % Determine number of timepoints in concatenated data
                    ntp2 = size(tss{cts}.ts{csub},2);                                                                                                            % Determine number of timepoints in current input
                    for cROI = 1:nROIs                                                                                                                           % For each ROI
                        tscat.ts{pm,1}(cROI,(ntp1+1):(ntp1+ntp2)) = tss{cts}.ts{csub}(cROI,:)-mean(tss{cts}.ts{csub}(cROI,:));                                   % Mean center each timeseries before concatenation
                    end
                case 2
                    ntp1 = size(tscat.ts{pm,1},2);                                                                                                               % Determine number of timepoints in concatenated data
                    ntp2 = size(tss{cts}.ts{csub},2);                                                                                                            % Determine number of timepoints in current input
                    for cROI = 1:nROIs                                                                                                                           % For each ROI
                        tscat.ts{pm,1}(cROI,(ntp1+1):(ntp1+ntp2)) = (tss{cts}.ts{csub}(cROI,:)-mean(tss{cts}.ts{csub}(cROI,:)))./std(tss{cts}.ts{csub}(cROI,:)); % z-score each timeseries before concatenation
                    end
            end
        end
    end
    
    [ds,dia] = setdiff(tss{cts}.subs,tscat.subs);                                                                                                          % Find participants present only in the current set
    for csub = 1:length(ds)                                                                                                                                % For each participant
        tscat.subs = [tscat.subs;ds{csub}];                                                                                                                % Set current identifier
        switch flag
            case 0
                tscat.ts = [tscat.ts;tss{cts}.ts{dia(csub)}];                                                                                              % Perform no manipulation
            case 1
                cind = length(tscat.subs);                                                                                                                 % Determine number of participants in currently concatenated data
                for cROI = 1:nROIs                                                                                                                         % For each ROI
                    tscat.ts{cind,1}(cROI,:) = tss{cts}.ts{dia(csub)}(cROI,:)-mean(tss{cts}.ts{dia(csub)}(cROI,:));                                        % Mean center each timeseries before concatenation
                end
            case 2
                cind = length(tscat.subs);                                                                                                                 % Determine number of participants in currently concatenated data
                for cROI = 1:nROIs                                                                                                                         % For each ROI
                    tscat.ts{cind,1}(cROI,:) = (tss{cts}.ts{dia(csub)}(cROI,:)-mean(tss{cts}.ts{dia(csub)}(cROI,:)))./std(tss{cts}.ts{dia(csub)}(cROI,:)); % z-score each timeseries before concatenation
                end
        end
    end
end