corr_measures_210802.m

%% ATK 210802 for loading from HD
% Calculate pairwise corr measures 
% Add mutual information

masterPath = "/Volumes/Aaron's PPC/ppc/2P_data/";
mouse = 'LD187';

%mySessions = {'LD187_141216','LD187_141215','LD187_141214','LD187_141213','LD187_141212',...
%    'LD187_141211','LD187_141210','LD187_141209','LD187_141208','LD187_141207','LD187_141206'};

mySessions = {'LD187_141216','LD187_141215','LD187_141214','LD187_141213'};
for i = 1:length(mySessions)
    session = cell2mat(mySessions(i));
    disp(session)

    % Load 2p data 
    output_dir_deconv = fullfile(masterPath, 'code_workspace',mouse,'deconvData');   
    sp = load(fullfile(output_dir_deconv,session),'c','b','c1','g','sn','sp').sp;
    syncedDataDir = fullfile(masterPath, 'code_workspace',mouse,'syncedData');   
    trialAlignedData = open(fullfile(syncedDataDir,[session '.mat'])).trialAlignedData;
    
    %% Calc pairwise mutual information
    % requires https://github.com/nmtimme/Neuroscience-Information-Theory-Toolbox
    %{
    
    % Format data for info theory toolbox and binarize activity
    % exclude first trial which has nans in ITI
    %allRaster = cell(1); allRaster{1} = sp>0; 
    lrRaster = cell(1); lrRaster{1} = trialAlignedData.Ca_concat(:,77:end) > 0;
    %trialMeanRaster = cell(1); trialMeanRaster{1} = trialAlignedData.Ca_trialMean_concat > 0;
    %trialResidualRaster = cell(1); trialResidualRaster{1} = trialAlignedData.Ca_residual_concat(:,77:end) > 0;
    
    % Note - trialResidualRaster is problematic because it is not always positive!

    tic
    numROIs = trialAlignedData.numCells;
    %pairMI_all = nan(numROIs, numROIs);
    pairMI_lr = nan(numROIs, numROIs);
    %pairMI_trialMean = nan(numROIs, numROIs);
    %pairMI_trialResidual = nan(numROIs, numROIs);
    parfor id1 = 1:numROIs
    %for id1 = 1:numROIs
        for id2 = 1:numROIs
            if id2 <= id1 % only calculate bottom triangle
                VariableIDs = {1,id1,1;1,id2,1}; % only relative timing important here
                %pairMI_all(id1,id2) = instinfo(allRaster, Method, VariableIDs);
                pairMI_lr(id1,id2) =  instinfo(lrRaster, 'PairMI', VariableIDs);
                %%%pairMI_trialMean(id1,id2) =  instinfo(trialMeanRaster, Method, VariableIDs);
                %pairMI_trialResidual(id1,id2) =  instinfo(trialResidualRaster, 'PairMI', VariableIDs);
            end
        end
        disp([num2str(id1) ' done']);
    end
    toc     
    %%
    pairMI_lr = tril(pairMI_lr,-1)' + tril(pairMI_lr);
    %pairMI_trialResidual = tril(pairMI_trialResidual,-1)' + tril(pairMI_trialResidual);
    
    %pairMI = struct;
    %pairMI.all = pairMI_all;
    %pairMI.lr = pairMI_lr;
    %pairMI.trialMean = pairMI_trialMean;
    %pairMI.trialResidual = pairMI_trialResidual;
    %}
    
    %% Calculate pairwise correlation matrices
    trialTypes = {'newL_trials','bR_trials','wL_trials','newR_trials'};
    pearsonCorr = struct;
    pearsonCorr.corr_all = corrcoef(sp');
    pearsonCorr.corr_Ca_lr = corrcoef(trialAlignedData.Ca_concat');
    pearsonCorr.corr_Ca_trialMean = corrcoef(trialAlignedData.Ca_trialMean_concat');
    pearsonCorr.corr_Ca_residual = corrcoef(trialAlignedData.Ca_residual_concat');
    
    % New metrics 220802
    % Trial mean averaged over b and w cues
    pearsonCorr.corr_trialMean_all =  corrcoef(trialAlignedData.trialMean_all');
    % Residual (activity minus trial mean above
    pearsonCorr.corr_trialResidual_bwavg =  corrcoef(trialAlignedData.trialResidual_bwavg_concat');
    % Trial mean diff bw b and w
    pearsonCorr.corr_bw_diff =  corrcoef(trialAlignedData.bw_diff');
    % Averaged over time, leave bw trial structure
    pearsonCorr.corr_timeMean =  corrcoef(trialAlignedData.timeMean');
    
    % New metrics 230127
    bin = trialAlignedData.Ca_concat>0; % use only timepoints during L/R trials
    pearsonCorr.corr_lr_bin = corrcoef(bin');
    dots = bin*bin';
    counts = sum(bin'); % total time points for each neuron with non-zero act
    counts_sum = bsxfun(@plus,counts,counts'); % avg of counts for pre and post neuron
    frac_coincident = dots./(counts_sum/2); % coincident counts / avg pre&post counts
    pearsonCorr.frac_coincident = frac_coincident;
    
    % Add also pair MI as comparison
    %pearsonCorr.pairMI = pairMI_lr;
    %% Save data as .mat
    output_dir =  fullfile(masterPath, 'code_workspace',mouse,'corrMetrics');   
    %save(fullfile(output_dir,session),'pearsonCorr', 'pairMI','-v7.3');
    save(fullfile(output_dir,session),'pearsonCorr','-v7.3');
end