CorrelationHeatmaps2022.m

%[OutputAverageMatrix, OutputNeuronOrder, OutputStructure]= CorrelationHeatmaps2022(MainStruct,InputType,FunctionalMeasure,Sorting,NeuronOrder)

% Inputs:
%
%       MainStruct, : This is the dataset structure where all the datasets
%       from a single condition (e.g. WT) are available. Compatible with
%       the structures generated by merging_head_and_tail.m 
%
%       InputType, : 'derivs' or 'traces',
%
%       FunctionalMeasure,   : 'Correlation', 'CrossCorrelation' or    'Covariance'
%    
%       Sorting,    : true or false, %if true the matrix will be sorted according
%       to hierarchical clustering, otherwise it will use the pre-sorted derivs order.
%
%       NeuronOrder, : List of the neurons (and the exact order, if sorting
%       is false) where the function searches within the datasets and
%       generates the average correlation matrix.
%
% Outputs:
%
%       OutputAverageMatrix: Average Correlation Matrix
%
%       OutputNeuronOrder: Neuron IDs and Order of the final Correlation Matrix
%        %if sorting is true it will give out the newly sorted version
%        %otherwise it is the same with the input neuron order
%
%       OutputStructure: This strcuture contains correlation matrices of
%       individual datasets, average correlation matrix across datasets and
%       the exact traces from each dataset from which subsequent
%       correlation matrices are calculated.

%example usage:
%[OutputAvMatr, OutputNeurOrd, OutputStruct]= CorrelationHeatmaps2022(Uzel_WT,'derivs','Correlation',0,NeuronOrder66);



%%
function [OutputAverageMatrix, OutputNeuronOrder, OutputStructure]= CorrelationHeatmaps2022(MainStruct,InputType,FunctionalMeasure,Sorting,NeuronOrder)

DatasetNumber=size(MainStruct,2); %this determines the number of datasets are available in the input structure

for n=1:DatasetNumber; % loop that will go through number of datasets within the input structure
    
    Number=num2str(n); %will be needed to label individual correlation matrices within the OutputStructure

    DatasetNeuronIDs=MainStruct(n).IDs; %taking the array with IDs from a particular dataset

    
    IDx=zeros(1,length(NeuronOrder)); %finding the indices of the neurons (given in the input NeuronOrder)
    for i=1:length(NeuronOrder); 
        for f=1:length(DatasetNeuronIDs);
            if strcmp(DatasetNeuronIDs{f}, NeuronOrder{i});
                IDx(i)=f;
            end
        end
    end
    
    
    Name1=strcat('InputTraces',Number); %IDed traces, dims: neuron number x frames
    Name2=strcat('FunctionalMatrix',Number);
    
        
    OutputStructure.(Name1)=zeros(size(MainStruct(n).traces,1),size(NeuronOrder,2));
    
    for j=1:length(NeuronOrder);
        if IDx(j)>0; %this means if the searched neuron is identified in the dataset, thus have a number other than 0
            
            switch InputType %taking the neuron trace either in original form or it's time-derivative
                case 'derivs'                                         
                    OutputStructure.(Name1)(:,j)=MainStruct(n).derivatives(:,IDx(j));
                case 'traces' 
                    OutputStructure.(Name1)(:,j)=MainStruct(n).traces(:,IDx(j));
            end
        else
            OutputStructure.(Name1)(:,j)=NaN; %if the searched neuron is not identified in the dataset, it will have NaNs instead.
        end
    end
    
    switch FunctionalMeasure %here the functional measure of choice wil be calculated from template structures
        case 'Correlation';
            
            [OutputStructure.(Name2)]=corrcoef(OutputStructure.(Name1));
            
            OutputStructure.preAVG(1:length(NeuronOrder),1:length(NeuronOrder),n)=OutputStructure.(Name2); %saving it as correlation matrix of the individual dataset 
            
        case 'Covariance';
            
            OutputStructure.(Name2)=cov(OutputStructure.(Name1));
            OutputStructure.preAVG(1:length(NeuronOrder),1:length(NeuronOrder),n)=OutputStructure.(Name2);

        case 'CrossCorrelation';
            
            tensecond=round((MainStruct(n).fps)*10); %finding the frame equivalent of ten seconds within the dataset
           
            
            Imatrix=nan(length(NeuronOrder)); %the matrix where the time lag is saved in case needed for further analysis
            for ki=1:length(NeuronOrder);
                for kj=1:length(NeuronOrder);
                    
                    if isnan(OutputStructure.(Name1)(1,ki)) || isnan(OutputStructure.(Name1)(1,kj))
                                 Xcorrmatr_temp(ki,kj)=nan;
                    else
                        

                   if max(xcorr(OutputStructure.(Name1)(:,ki),OutputStructure.(Name1)(:,kj),tensecond,'coeff'))==max(abs(xcorr(OutputStructure.(Name1)(:,ki),OutputStructure.(Name1)(:,kj),tensecond,'coeff')));
                      
                        [Xcorrmatr_temp(ki,kj),Itemp]=max(xcorr(OutputStructure.(Name1)(:,ki),OutputStructure.(Name1)(:,kj),tensecond,'coeff'));
                 
                        Imatrix(ki,kj)=abs(Itemp-length(xcorr(OutputStructure.(Name1)(:,ki),OutputStructure.(Name1)(:,kj),tensecond,'coeff'))/2)/MainStruct.(DatasetNumber).fps;
                   else
                       [Xcorrmatr_temp(ki,kj),Itemp]=min(xcorr(OutputStructure.(Name1)(:,ki),OutputStructure.(Name1)(:,kj),tensecond,'coeff'));
                        Imatrix(ki,kj)=abs(Itemp-length(xcorr(OutputStructure.(Name1)(:,ki),OutputStructure.(Name1)(:,kj),tensecond,'coeff'))/2)/MainStruct.(DatasetNumber).fps;
                        
                   end
                    
                    end
                    
                end
            end
            OutputStructure.Imatrs(1:length(NeuronOrder),1:length(NeuronOrder),n)=Imatrix;
            OutputStructure.(Name2)=Xcorrmatr_temp;
            OutputStructure.preAVG(1:length(NeuronOrder),1:length(NeuronOrder),n)=OutputStructure.(Name2);
    end
    clear IDx ID
    clearvars -except Sorting  OutputStructure InputType FunctionalMeasure MainStruct NeuronOrder  
end

%%
OutputStructure.AverageMatr=nanmean(OutputStructure.preAVG,3); %calculating the average matrix across all datasets

if Sorting==true
    
    LinkageMatrix=OutputStructure.AverageMatr; %generating the variable LinkageMAtrix for sorting purposes
    LinkageMatrix(isnan(LinkageMatrix)) = 0;
    NumTracks=size(LinkageMatrix,1);
    
    [~,~,outperm]=dendrogram(linkage(LinkageMatrix),NumTracks); %outperm gives out the sorted order.
    close all
    %
    
 NeuronOrder_resorted=NeuronOrder(outperm); %rearranging NeuronOrder
 

        
    for i=1:length(NeuronOrder);
        for j=1:length(NeuronOrder);
            
            newindices=[outperm(i),outperm(j)];
            ReSorted_matrix_HIS(i,j)=OutputStructure.AverageMatr(newindices(1),newindices(2));
            
        end
    end
    
    OutputStructure.AverageMatr=ReSorted_matrix_HIS;
    
    OutputNeuronOrder=NeuronOrder_resorted;
end


OutputAverageMatrix=OutputStructure.AverageMatr;
OutputNeuronOrder=NeuronOrder;

end