Merge branch 'tcca_glm_v3'

% 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, rcMap, glmSolveMethod, idxBasis, paramsBasis, rhoSD_ssThresh, flagNuisanceRMethod, driftOrder, flagMotionCorrect )

%

% UI NAME:

% GLM_HRF_Drift_SS

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

% tIncAuto - a vector (#time points x 1) indicating which data time points

%            are motion (=0) or not (=1)

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

% trange - defines the range for the block average [tPre tPost]

% glmSolveMethod - this specifies the GLM solution method to use

%            1 - use ordinary least squares (Ye et al (2009). NeuroImage, 44(2), 428?447.)

%          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

%

%

%

% 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, rcMap, trange, glmSolveMethod, idxBasis, paramsBasis, rhoSD_ssThresh, flagNuisanceRMethod, driftOrder, flagMotionCorrect)

%

%

% 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

%

%

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, rcMap, trange, glmSolveMethod, idxBasis, paramsBasis, rhoSD_ssThresh, flagNuisanceRMethod, driftOrder, flagMotionCorrect)

% Init output

data_yavg     = DataClass().empty();

    end

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

    if isempty(lstSS)

    if isempty(lstSS) || (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 ischar(rcMap{iBlk}) % use all regressors for all channels (only one group)

                    mlSSlst = 1;

                elseif iscell(rcMap{iBlk}) % 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 {0,1,2} % short separation

            for iConc=1:2

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

            end

        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

    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 ischar(rcMap{iBlk}) % no channel map: use all tCCA regressors for one group of all channels

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

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

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

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

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

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

                end

            end

            if ~isempty(lstML)

% SYNTAX:

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

% [Aaux, rcMap] = hmrR_tCCA(data, aux, probe, runIdx, flagtCCA, tCCAparams, tCCAaux_inx, rhoSD_ssThresh, runIdxResting, tResting)

%

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

% CITATION:

% Canonical Correlation Anlaysis described in detail in von Lhmann et al. (2020), NeuroImage. 

% If you use this method, please use the following CITATION:

% von Lhmann, Alexander, et al. "Improved physiological noise regression

% in fNIRS: A multimodal extension of the General Linear Model using temporally

% embedded Canonical Correlation Analysis." NeuroImage 208 (2020): 116472.

% probe - SNIRF probe type containing source/detector geometry data (See SNIRF Spec for more details)

% runIdx - the index of the run in a multi-run session

% flagtCCA - turns the function on / off

% flagICRegressors - selects regressor generation strategy:

%            0/false - common set of regressors for all fNIRS channels (default)

%            1/true - one individual regressor per fNIRS channel (true)

% tCCAparams - These are the parameters for tCCA function

%            1 - timelag

%            2 - sts: step size

%            3 - ctr: if <1  correlation threshold, if =>1 number of regressors.

%            1 - 1 timelag [s]

%            2 - 2 step size [samples]

%            3 - 3 ctr if <1  correlation threshold, if =>1 number of regressors.

%                   Redundant if flagICRegressors = 1;

% tCCAaux_inx -  Indices of the aux channels to be used to generate regressors

% rhoSD_ssThresh - max distance for a short separation measurement. Set =0

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

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

% runIdxResting - resting state run index

% tResting - start/stop time [s] to use from resting data for tCCA training

%

% OUTPUTS:

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

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

%           aux regressor indices for individual regressor-channel mapping.

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

% rcMap - Currently always empty or "all". 

%           Under development: Will also provide an array of cells (1 x #fNIRS channels) 

%           containing aux regressor indices for individual regressor-channel mapping.

%           Only relevant when flagICRegressors = 1.

%

% USAGE OPTIONS:

% hmrR_tCCA_Concentration_Data: [Aaux, rcMap] = hmrR_tCCA(dc, aux, probe, iRun, flagtCCA, flagICRegressors, tCCAparams, tCCAaux_inx, rhoSD_ssThresh, runIdxResting)

% hmrR_tCCA_OD_Data: [Aaux, rcMap] = hmrR_tCCA(dod, aux, probe, iRun, flagtCCA, flagICRegressors, tCCAparams, tCCAaux_inx, rhoSD_ssThresh, runIdxResting)

% hmrR_tCCA_Concentration_Data: [Aaux, rcMap] = hmrR_tCCA(dc, aux, probe, iRun, flagtCCA, tCCAparams, tCCAaux_inx, rhoSD_ssThresh, runIdxResting, tResting)

% hmrR_tCCA_OD_Data: [Aaux, rcMap] = hmrR_tCCA(dod, aux, probe, iRun, flagtCCA, tCCAparams, tCCAaux_inx, rhoSD_ssThresh, runIdxResting, tResting)

%

%

% PARAMETERS:

% flagtCCA: 1

% flagICRegressors: 0

% tCCAparams: [3 2 0.3]

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

% rhoSD_ssThresh: 15.0

% runIdxResting: 1

% tResting: [30 210]

%

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

function [Aaux, rcMap] = hmrR_tCCA(data, aux, probe, runIdx, flagtCCA, tCCAparams, tCCAaux_inx, rhoSD_ssThresh, runIdxResting, tResting)

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

%%

%% flags and tCCA settings

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

% flags.shrink = true; % perform shrinkage in the CCA

flags.shrink = false;  % JD: set to false because rtcca generate error, "Undefined function or variable 'cshrink'".

% flagICRegressors - selects regressor generation strategy. (0/false) chooses a common set

% of regressors for all fNIRS channels. (1/true) generates one individual

% regressor per fNIRS channel. This feature is currently under development,

% please do not use it.

flagICRegressors = 0;

% tCCA parameters

param.tau = tCCAparams(2); %stepsize for embedding in samples (tune to sample frequency!)

timelag = tCCAparams(1);

param.ct = 0;   % correlation threshold for rtcca function, set here to 0 (will be applied later)

% correlation used outside of the rtcc function

% correlation used outside of the rtcca function

ctr = tCCAparams(3);

%checksum for tCCA filter file

chksm = int16(7*flagICRegressors + sum(11*tCCAparams(:)) + sum(13*tCCAaux_inx(:)) + 17*rhoSD_ssThresh + 19*runIdxResting + sum(23*tResting(:)));

if flagtCCA

    SrcPos  = probe.GetSrcPos();

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

        %% DO THE TCCA WORK

        filterFilename = sprintf('./tCCAfilter_%d.txt', iBlk);

        tCCAexists = isfile(sprintf('./tCCAfilter_%d.txt', iBlk));

        filterFilename = sprintf('./tCCAfilter_%d_%d.txt', iBlk, chksm);

        tCCAexists = isfile(sprintf('./tCCAfilter_%d_%d.txt', iBlk, chksm));

        isTrainingRun = runIdxResting == runIdx;

        if ~tCCAexists && isTrainingRun

            dotCCA = 'train';

        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

                %       Columns of the tCCA filter matrix correspond to common regressors

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

                %       number of embeddings

                if flagICRegressors %% learn channel specific filters and regressors

                    for cc = 1:size(d_long,2)

                        %% perform tCCA with shrinkage for each fNIRS channel

                        [REG,  ADD] = rtcca_dummy(d_long(:,cc),AUX,param,flags);

                        %% save individual regressor and channel map

                        % return the reduced number of available regressors

                        Aaux{iBlk}(:,cc) = REG(:,1);

                        % assemble tCCA filter for channel specific regressors from individual reduced mapping matrices Av

                        tCCAfilter(:,cc) = ADD.Av(:,1);

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

                        rcMap{iBlk}{cc} = cc;

                %       number of embeddings. If flagICRegressors = 1 (currently N/A), 

                %       each column corresponds to an indiviudal channel regressor. 

                % cut data to selected time window before training

                % warning if resting segment is shorter than 1min

                if diff(tResting)<60

                    msgbox('WARNING: tCCA training with less than 60s of data is not recommended!')

                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)

                cstrt = find(t>=tResting(1));

                cstp = find(t<=tResting(2));

                cIdx = [cstrt(1):cstp(end)];

                % cut

                d_long = d_long(cIdx,:);

                AUX = AUX(cIdx,:);

                if ~flagICRegressors %% learn common set of regressors for all channels (default)

                    %% perform tCCA with shrinkage for all fNIRS channels

                    [REG,  ADD_trn] = rtcca(d_long,AUX,param,flags);                 %% save and output tCCA filter.

                    %reduce filter matrix with the help of correlation threshold or max number of regressors

                    if ctr < 1

                        % use only auxiliary tcca components that have correlation > ct

                        compindex=find(ADD_trn.ccac>param.ct);

                        %% JD had a comment here as Undefined variable "ADD_trn" or class "ADD_trn.ccac"

                        compindex=find(ADD_trn.ccac > ctr);

                        % throw a warning that overfitting might occur if

                        % all regressor's correlations are > ctr

                        if numel(compindex) == numel(ADD_trn.ccac)

                            msgbox('WARNING: All regressors have a canonical correlation larger than your threshold. To avoid overfitting, consider setting the number of regressors to a fixed number ctr = N.')

                        end

                    else

                        % use only the first ctr auxiliary tcca components (fixed number of

                        % regressors = ctr)

                    end

                    % return the reduced number of available regressors

                    Aaux{iBlk} = REG(:,compindex);

                    Aaux = REG(:,compindex);

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

                    tCCAfilter = ADD_trn.Av(:,compindex);

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

                    rcMap{iBlk} = 'all';

                else %% learn channel specific filters and regressors. UNDER DEVELOPMENT - DO NOT USE

                    for cc = 1:size(d_long,2)/2

                        %% perform tCCA with shrinkage for each fNIRS channel (always HbO and HbR together)

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

                        %% save individual regressor and channel map

                        % return the reduced number of available regressors

                        Aaux{iBlk}(:,cc) = REG(:,1);

                        % assemble tCCA filter for channel specific regressors from individual reduced mapping matrices Av

                        tCCAfilter(:,cc) = ADD.Av(:,1);

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

                        rcMap{iBlk}{cc} = cc;

                    end

                    % reshape (HbO first row, HbR second row)

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

                end

                Aaux = []; % we do not want to use regressors for the resting run (available only for investigative purposes)

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

                fprintf('hmrR_tCCA: run idx = %d. Generated and Saved tCCAfilter\n', runIdx)

                save(filterFilename, '-ascii', 'tCCAfilter');

                fprintf('Canonical correlation coefficients of all trained regressors:')

                ADD_trn.ccac(compindex)

            case 'apply'

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

                % Load the filter for the iBlk data block

            case 'skip'

                Aaux = [];

                rcMap = [];

                %% put a user warning

                if runIdx == runIdxResting-1

                    msgbox('tCCA raining (resting run) is not the first run. Other runs skipped. Please re-run the session for complete results.')

                elseif runIdx > runIdxResting

                    msgbox('no tCCA filter trained. Please run training resting run or whole session first.')

                end

        end

    end