updated tCCA and added tCCA support in GLM function (first try)

% SYNTAX:

% [yavg, yavgstd, tHRF, nTrials, ynew, yresid, ysum2, beta, R] = hmrR_GLM(data, stim, probe, mlActAuto, Aaux, tIncAuto, trange, glmSolveMethod, idxBasis, paramsBasis, rhoSD_ssThresh, flagSSmethod, driftOrder, flagMotionCorrect )

% [yavg, yavgstd, tHRF, nTrials, ynew, yresid, ysum2, beta, R] = hmrR_GLM(data, stim, probe, mlActAuto, Aaux, tIncAuto, trange, glmSolveMethod, idxBasis, paramsBasis, rhoSD_ssThresh, flagNuisanceRMethod, driftOrder, flagMotionCorrect, rcMap )

%

% UI NAME:

% GLM_HRF_Drift_SS

%          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.

% flagSSmethod - 0 if short separation regression is performed with the nearest

% flagNuisanceRMethod - 0 if short separation regression is performed with the nearest

%               short separation channel.

%            1 if performed with the short separation channel with the

%               greatest correlation.

%            2 if performed with average of all short separation channels.

%            3 uses tCCA regressors for nuisance regression, in Aaux,

%            mapped by rcMap, provided by hmr_tCCA()

% driftOrder - Polynomial drift correction of this order

% flagMotionCorrect - set to 1 to baseline correct between motion epochs indicated in tIncAuto, otherwise set to 0

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

%           indices for individual regressor-channel mapping. Currently

%           only relevant when flagNuisanceRMethod = 3 (tCCA regressors).

%

%

% OUTPUTS:

%

%

% USAGE OPTIONS:

% GLM_HRF_Drift_SS_Concentration: [dcAvg, dcAvgStd, nTrials, dcNew, dcResid, dcSum2, beta, R] = hmrR_GLM(dc, stim, probe, mlActAuto, Aaux, tIncAuto, trange, glmSolveMethod, idxBasis, paramsBasis, rhoSD_ssThresh, flagSSmethod, driftOrder, flagMotionCorrect )

% GLM_HRF_Drift_SS_Concentration: [dcAvg, dcAvgStd, nTrials, dcNew, dcResid, dcSum2, beta, R] = hmrR_GLM(dc, stim, probe, mlActAuto, Aaux, tIncAuto, trange, glmSolveMethod, idxBasis, paramsBasis, rhoSD_ssThresh, flagNuisanceRMethod, driftOrder, flagMotionCorrect, rcMap )

%

%

% PARAMETERS:

% idxBasis: 2

% paramsBasis: [0.1 3.0 10.0 1.8 3.0 10.0], maxnum: 6

% rhoSD_ssThresh: 15.0

% flagSSmethod: 1

% flagNuisanceRMethod: 1

% driftOrder: 3

% flagMotionCorrect: 0

% rcMap: []

%

%

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

    hmrR_GLM(data_y, stim, probe, mlActAuto, Aaux, tIncAuto, trange, glmSolveMethod, idxBasis, paramsBasis, rhoSD_ssThresh, flagSSmethod, driftOrder, flagMotionCorrect )

    hmrR_GLM(data_y, stim, probe, mlActAuto, Aaux, tIncAuto, trange, glmSolveMethod, idxBasis, paramsBasis, rhoSD_ssThresh, flagNuisanceRMethod, driftOrder, flagMotionCorrect, rcMap )

% Init output

data_yavg     = DataClass().empty();

    end

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

    if isempty(lstSS)

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

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

        mlSSlst = 0;

    else

        if flagSSmethod==0  % use nearest SS

        switch flagNuisanceRMethod

            case 0  % use nearest SS

                for iML = 1:length(lst)

                    rho = sum((ones(length(lstSS),1)*posM(iML,:) - posM(lstSS,:)).^2,2).^0.5;

                    [foo,ii] = min(rho);

                    iNearestSS(iML) = lstSS(ii);

                end

                mlSSlst = unique(iNearestSS);       

        elseif flagSSmethod==1 % use SS with highest correlation

            case 1 % use SS with highest correlation    

                % HbO

                dc = squeeze(y(:,1,:));

                dc = (dc-ones(length(dc),1)*mean(dc,1))./(ones(length(dc),1)*std(dc,[],1));

                cc(:,:,2) = dc'*dc / length(dc);

                clear dc           

                % find short separation channel with highest correlation

                for iML = 1:size(cc,1)

                    % HbO

                    [foo,ii] = max(cc(iML,lstSS,2));

                    iNearestSS(iML,2) = lstSS(ii);

                end 

        elseif flagSSmethod==2 % use average of all active SS as regressor

            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)

                    mlSSlst = 1;

                else % use channel regressor map

                    mlSSlst = 1:size(rcMap,2); 

                end

        end

    end

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

    % Expand design matrix with Aaux

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

    nAux = size(Aaux,2);

    if flagNuisanceRMethod == 3

        nAux = 0;

    end

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

    % Expand design matrix for Motion Correction

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

    % Final design matrix

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

    switch flagNuisanceRMethod

        case {1,2,3} % 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)

            for iConc=1:2

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

            end

    end

    nCh = size(y,3);

    foo     = zeros(nB*nCond+driftOrder+1+nAux+nMC,nCh,2); % 4 extra for 3rd order drift + nAux

    for conc=1:2 %only HbO and HbR

        if flagSSmethod==1 & ~isempty(lstSS) %rhoSD_ssThresh>0

        if flagNuisanceRMethod==1 & ~isempty(lstSS) %rhoSD_ssThresh>0

            mlSSlst = unique(iNearestSS(:,conc));

        end

        % loop over short separation groups

        % loop over short separation / nuisance regressor groups

        for iSS = 1:length(mlSSlst)

            lstMLtmp = 1:size(ml,1);

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

                % lstML = 1:size(y,3);

                At = A(:,:,conc);

            elseif flagSSmethod==0

            elseif flagNuisanceRMethod==0

                lstML = find(iNearestSS(:)==mlSSlst(iSS) & mlAct(lstMLtmp)==1);

                % lstML = find(iNearestSS==mlSSlst(iSS));

                Ass = y(:,conc,mlSSlst(iSS));

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

            elseif flagSSmethod==1

            elseif flagNuisanceRMethod==1

                lstML = find(iNearestSS(:,conc)==mlSSlst(iSS) & mlAct(lstMLtmp)==1);

                % lstML = find(iNearestSS(:,conc)==mlSSlst(iSS));

                Ass = y(:,conc,mlSSlst(iSS));

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

            elseif flagSSmethod==2

            elseif flagNuisanceRMethod==2

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

                % lstML = 1:size(y,3);

                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

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

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

                else % 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];

                end

            end

            if ~isempty(lstML)

% 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

%

%          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 - filter matrix that was learned from resting state data.

%           Will be freshly trained if input was empty. 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

% 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 column corresponds to an indiviudal channel

%           regressor. Otherwise, columns correspond to common regressors

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

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

% 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