Commit 95396211 authored by mayucel's avatar mayucel
Browse files

incorporating Alex and Jay code

needs to be triple checked!
parent 529848ce
Loading
Loading
Loading
Loading
+106 −92
Original line number Diff line number Diff line
@@ -2,7 +2,7 @@
% [Aaux, rcMap] = hmrR_tCCA(data, aux, probe, runIdx, flagtCCA, flagICRegressors, tCCAparams, tCCAaux_inx, rhoSD_ssThresh, runIdxResting)
%
% UI NAME:
% User_Friendly_Name_For_hmrR_tCCA
% tCCA regressor generator
%
% DESCRIPTION:
% This script generates regressors using the regularized temporally embedded
@@ -139,16 +139,32 @@ if flagtCCA
        % zscore AUX signals
        AUX = zscore(AUX);
        
        param.NumOfEmb = ceil(timelag*fq / tCCAparams(2));

        %% DO THE TCCA WORK
        filterFilename = sprintf('./tCCAfilter_%d.txt', iBlk);
        tCCAexists = isfile(sprintf('./tCCAfilter_%d.txt', iBlk));
        isTrainingRun = runIdxResting == runIdx;
        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 / param.tau);
        if runIdx == runIdxResting
            %% if the tCCAfilter variable is empty, learn and put it out (this is the training/resting run)
                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);    %%%%%%% JD: rtcca generates error, "Undefined function or variable 'cshrink'". %%%%%%% 
                        [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);
@@ -157,48 +173,45 @@ if flagtCCA
                        % save channel-regressor map. First half HbO, second half HbR
                        rcMap{iBlk}{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_dummy(d_long,AUX,param,flags);    %%%%%%% JD: rtcca generates error, "Undefined function or variable 'cshrink'". %%%%%%% 
                %% save and output tCCA filter.
                    [REG,  ADD] = rtcca_dummy(d_long,AUX,param,flags);                 %% save and output tCCA filter.
                    %reduce filter matrix with the help of correlation threshold or max number of regressors
                
                %%%%%%% JD: Undefined variable "ADD_trn" or class "ADD_trn.ccac". %%%%%%% 
                    if ctr < 1
                        % use only auxiliary tcca components that have correlation > ct
                    % compindex=find(ADD_trn.ccac>param.ct);
                    compindex=1:10;
                        compindex=find(ADD_trn.ccac>param.ct);
                        %% JD had a comment here as Undefined variable "ADD_trn" or class "ADD_trn.ccac"
                    else
                        % use only the first ctr auxiliary tcca components (fixed number of
                        % regressors = ctr)
                    % if numel(ADD_trn.ccac) >= ctr
                    %     compindex = 1:ctr;
                    % else % not as many components available as regressors requested, provide all
                    %     compindex = 1:numel(ADD_trn.ccac);
                    % end                    
                    compindex=1:10;
                        if numel(ADD_trn.ccac) >= ctr
                            compindex = 1:ctr;
                        else % not as many components available as regressors requested, provide all
                            compindex = 1:numel(ADD_trn.ccac);
                        end
                        
                    end
                    % return the reduced number of available regressors
                    Aaux{iBlk} = REG(:,compindex);
                    % 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{iBlk} = [];
                    rcMap{iBlk} = 'all';
                end
                %% save tCCAfilter matrix that was learned from resting state data.
                fprintf('hmrR_tCCA: run idx = %d. Generated and Saved tCCAfilter\n', runIdx)
                print_filter(tCCAfilter);
                save(filterFilename, '-ascii', 'tCCAfilter');
                
        elseif exist(filterFilename,'file')
            
            %% if the tCCAfilter variable exists, apply the filtering and generate the tCCA regressors
            fprintf('hmrR_tCCA: run idx = %d. Loading and Using tCCAfilter\n', runIdx)
            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
                fprintf('hmrR_tCCA: run idx = %d. Loading and Using tCCAfilter\n', runIdx)
                tCCAfilter = load(filterFilename,'-ascii');
                print_filter(tCCAfilter);
                
            
                % Temporal embedding and zscoring of auxiliary data
                aux_sigs = AUX;
                aux_emb = aux_sigs;
@@ -208,29 +221,30 @@ if flagtCCA
                    aux_emb = [aux_emb aux];
                end
                %zscore
            % aux_emb = zscore(aux_emb);
            aux_emb = 1;     % JD: Dummy value 
                aux_emb = zscore(aux_emb);
                % Apply tCCA filter (generate regressors)
            Aaux{iBlk} = aux_emb * tCCAfilter;
                Aaux = aux_emb * tCCAfilter;
                % 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
                        rcMap{iBlk}{cc} = cc;
                    end
                    % reshape (HbO first row, HbR second row)
                    rcMap = reshape(rcMap,[2,numel(rcMap)/2]);
                else
                rcMap{iBlk} = [];
                    rcMap{iBlk} = 'all';
                end
        
        else
            
            Aaux{iBlk} = [];
            rcMap{iBlk} = [];
            % tCCAfilter = []; if uncommented, retraining necessary whenever tcca is temporarily flagged out of processing stream
            
            case 'skip'
                Aaux = [];
                rcMap = [];
        end
        
    end
else
    Aaux = [];
    rcMap = [];
    % tCCAfilter = []; if uncommented, retraining necessary whenever tcca is temporarily flagged out of processing stream
end