minor debugs and moved functions to main UserFucntions folder for testing

% flagNuisanceRMethod: 1

% driftOrder: 3

% flagMotionCorrect: 0

% rcMap: []

% rcMap: 'all'

%

%

function [data_yavg, data_yavgstd, nTrials, data_ynew, data_yresid, data_ysum2, beta_blks, yR_blks] = ...

    end

    lstSS = lst(find(rhoSD<=rhoSD_ssThresh & mlAct(lst)==1));

    if isempty(lstSS) || (isempty(rcMap) && isempty(Aaux) && flagNuisanceRMethod == 3)

    if isempty(lstSS) || (isempty(Aaux) && flagNuisanceRMethod == 3)

        fprintf('There are no short separation channels in this probe ...performing regular deconvolution.\n');

        mlSSlst = 0;

    else

            case 2 % use average of all active SS as regressor

                mlSSlst = 1;

            case 3 % use tCCA regressors and channel map from hmrR_tCCA()

                if isempty(rcMap) % use all regressors for all channels (only one group)

                if ischar(rcMap) % use all regressors for all channels (only one group)

                    mlSSlst = 1;

                else % use channel regressor map

                elseif iscell(rcMap) % use channel regressor map

                    mlSSlst = 1:size(rcMap,2); 

                end

    % Final design matrix

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    switch flagNuisanceRMethod

        case {1,2,3} % short separation

        case {0,1,2} % short separation

            for iConc=1:2

                A(:,:,iConc)=[dA(:,:,iConc) xDrift Aaux Amotion];

            end

        case 4 % tCCA regressor design matrix without Aaux (will be put in in the loop below)

        case 3 % tCCA regressor design matrix without Aaux (will be put in in the loop below)

            for iConc=1:2

                A(:,:,iConc)=[dA(:,:,iConc) xDrift Amotion];

            end

                Ass = mean(y(:,conc,lstSS),3);

                At = [A(:,:,conc) Ass];

            elseif flagNuisanceRMethod==3

                if isempty(rcMap) % no channel map: use all tCCA regressors for one group of all channels

                if ischar(rcMap) % no channel map: use all tCCA regressors for one group of all channels

                    lstML = lstMLtmp(find(mlAct(lstMLtmp)==1));

                    At = [A(:,:,conc) Aaux]; 

                else % channel map: each single regressor corresponds to one channel (nCH groups)

                elseif iscell(rcMap) % channel map: each single regressor corresponds to one channel (nCH groups)

                    lstML = lstMLtmp(find(mlAct(rcMap{conc,iSS})==1));

                    Atcca = Aaux(:,rcMap{conc,iSS});

                    At = [A(:,:,conc) Atcca];

% SYNTAX:

% [Aaux, tCCAfilter, rcMap] = hmrR_tCCA(data, aux, probe, flagtCCA, flagICRegressors, tCCAparams, tCCAaux_inx, rhoSD_ssThresh, tCCAfilter)

% [Aaux, rcMap] = hmrR_tCCA(data, aux, probe, flagtCCA, flagICRegressors, tCCAparams, tCCAaux_inx, rhoSD_ssThresh, tCCAtrnRunIdx, runIndex)

% UI NAME:

% hmrR_tCCA

%

% DESCRIPTION:

% This script generates regressors using the regularized temporally embedded

% Canonical Correlation Anlaysis described in detail in von Lhmann, NI,

% 2020. Generated Aaux is to be passed as an input to GLM function.

% 2020. Generated Aaux and rcMap are to be passed as an input to GLM function.

% Reference code:  https://github.com/avolu/tCCA-GLM

%

% If you apply this method or want to learn more about it, please use the following

%          if you do not want to regress the short separation measurements.

%          Follows the static estimate procedure described in Gagnon et al (2011).

%          NeuroImage, 56(3), 1362?1371.

% tCCAfilter - cell array of the filter matrix/matrices that was/were learned from resting state data.

%           Will be freshly trained if input is empty. If flagICRegressors = 1, 

%           then the dimension of the cell array is (1 x #fNIRS channels), otherwise (1x1).

% tCCAtrnRunIdx - index of the run used for tCCA training (resting state

%          run), user selection

% runIndex - index of the current run (provided internally)

%

% OUTPUTS:

% Aaux - A matrix of auxilliary regressors (#time points x #Aux regressors)

% tCCAfilter - filter matrix that was learned from resting state data.

%           Will be freshly trained if input was empty. If flagICRegressors = 1, 

%           then each regressor 

% rcMap - An array of cells (1 x #fNIRS channels) containing

%           for all channels in descending order ranked by canonical correlation coefficient

% rcMap - An array of cells (2 (HbO+HbR) x #fNIRS channels) containing

%           aux regressor indices for individual regressor-channel mapping. 

%           Only relevant when flagICRegressors = 1, otherwise rcMap is empty.

%

%

% USAGE OPTIONS:

% [Aaux, tCCAfilter, rcMap] = hmrR_tCCA(data, aux, probe, flagtCCA, flagICRegressors, tCCAparams, tCCAaux_inx, rhoSD_ssThresh, tCCAfilter)

% [Aaux, rcMap] = hmrR_tCCA(data, aux, probe, flagtCCA, flagICRegressors, tCCAparams, tCCAaux_inx, rhoSD_ssThresh, tCCAtrnRunIdx, runIndex)

%

%

% PARAMETERS:

% tCCAparams: [3 2 0.3] 

% tCCAaux_inx: [1 2 3 4 5 6 7 8]

% rhoSD_ssThresh: 15.0

% tCCAfilter: []

% tCCAtrnRunIdx: 1

%% COMMENTS/THOUGHTS/QUESTIONS ALEX

% 2) Output canonical correlation coefficients as quality metric? 

% 3) Output fNIRS signal(s)/Aux regressors for visualization/quality control?

% 4) single channel (regressor) feature: regressor channel map currently assumes first half HbO, secon half HbR. CHECK HOW NEEDS TO BE ADAPTED FOR GLM?

% 5) Implement the variable low/bandpass filter coefficients from previous processing stream !!!

% 6) check/provide group session root directory

%%

function [Aaux, tCCAfilter, rcMap] = hmrR_tCCA(data, aux, probe, flagtCCA, flagICRegressors, tCCAparams, tCCAaux_inx, rhoSD_ssThresh, tCCAfilter)

function [Aaux, rcMap] = hmrR_tCCA(data, aux, probe, flagtCCA, flagICRegressors, tCCAparams, tCCAaux_inx, rhoSD_ssThresh, tCCAtrnRunIdx, runIndex)

%% flags and tCCA settings

flags.pcaf =  [0 0]; % no pca of X or AUX

    % only selected signals

    AUX = AUX(:,tCCAaux_inx);

    % lowpass filter AUX signals

    AUX = hmrR_BandpassFilt(AUX, fq, 0, 0.5);

%     AUX = hmrR_BandpassFilt(AUX, fq, 0, 0.5); % now in hmrR_bandpass

    % AUX signals + add ss signal if it exists

    if ~isempty(d_short)

        AUX = [AUX, d_short]; % this should be of the current run

    AUX = zscore(AUX);

    %% DO THE TCCA WORK

    if isempty(tCCAfilter)

        %% if the tCCAfilter variable is empty, learn and put it out (this is the training/resting run)

    % check whether filter has been trained, and exists in group root

    % directory, should be trained, or trial should be skipped

    tCCAexists = isfile('tCCAfilter.mat');

    isTrainingRun = tCCAtrnRunIdx == runIndex;

    if tCCAexists && isTrainingRun

        dotCCA = 'train';

    elseif tCCAexists && ~isTrainingRun

        dotCCA = 'apply';

    else

        dotCCA = 'skip';

    end

    switch dotCCA

        case 'train'

            %% if the tCCAfilter variable is not existing, learn and save it (this is the training/resting run)

            %       If flagICRegressors = 1,

            %       then each column corresponds to an indiviudal channel

            %       regressor. Otherwise, columns correspond to common regressors

            %       for all channels in descending order ranked by canonical correlation coefficient

            %       number of embeddings

            param.NumOfEmb = ceil(timelag*fq / sts);

            if flagICRegressors %% learn channel specific filters and regressors

                    % save channel-regressor map. First half HbO, second half HbR

                    rcMap{cc} = cc;

                end

                % reshape (HbO first row, HbR second row)

                rcMap = reshape(rcMap,[2,numel(rcMap)/2]);

            else %% learn common set of regressors for all channels (default)

                %% perform tCCA with shrinkage for all fNIRS channels

                [REG,  ADD] = rtcca(d_long,AUX,param,flags);

                % set channel-regressor map to empty (GLM will use all available regressors for all channels)

                rcMap = [];

            end

    else

        %% if the tCCAfilter variable exists, apply the filtering and generate the tCCA regressors

            %% save tCCAfilter matrix that was learned from resting state data.

            save('tCCAfilter', 'tCCAFilter');

        case 'apply'

            %% if the tCCAfilter variable exists, load it, apply the filtering and generate the tCCA regressors

            % load filter from group root directory

            load('tCCAfilter.mat')

            % Temporal embedding and zscoring of auxiliary data

            aux_sigs = AUX;

            aux_emb = aux_sigs;

            % if the individual channel regressor flag was set, indicate by providing rc mapping

            if flagICRegressors

                for cc = 1:size(d_long,2)

                % save channel-regressor map. First half HbO, second half HbR

                    % save channel-regressor map. First column HbO, second column HbR

                    rcMap{cc}=cc;

                end

                % reshape (HbO first row, HbR second row)

                rcMap = reshape(rcMap,[2,numel(rcMap)/2]);

            else

                rcMap = [];

            end

        case 'skip'

            Aaux = [];

            rcMap = [];

    end

else

    Aaux = [];

    rcMap = [];

    % tCCAfilter = []; if uncommented, retraining necessary whenever tcca is temporarily flagged out of processing stream

end

                % return reduced mapping matrix Av, this is the tCCA filter

                tCCAfilter = ADD.Av(:,compindex);

                % set channel-regressor map to empty (GLM will use all available regressors for all channels)

                rcMap = [];

                rcMap = 'all';

            end

            %% save tCCAfilter matrix that was learned from resting state data.

            save('tCCAfilter', 'tCCAFilter');

                % reshape (HbO first row, HbR second row)

                rcMap = reshape(rcMap,[2,numel(rcMap)/2]);

            else

                rcMap = [];

                rcMap = 'all';

            end

        case 'skip'

            Aaux = [];