%This is a script estimating time-delayed phase MI for a series of trials and provides an average 
%MI table for all trials of each condition, each timepoint, and each subject. The
%folder is expected to have a sub-directory with the files needed for the
%script to run called 'FilesNeeded', another one in which the output data will be saved called 'MusicPlast_WithinBetweenMI',
%the script it self, and the timeseries to be analyzed.
%Theoretically once the folder is the current directory of matlab, you may
%just run it with no further input needed. Just hit the green button and
%wait.... a lot!

%Add necessary toolbox GCMI to the path!! it is supposed to be in a sub-folder of the current directory called FilesNeeded
% addpath('./FilesNeeded/GCMI');
% create a list with all the files to be processed
filelist=dir('*_Hilbert.mat');

% Get lists of subjects, times and conditions to classify the files
filedetails = getdetails(filelist); %This calls a function placed at the end of the script that generates a table with the files' subject names, condition1, condition2.

%Load lists of files characteristics to be used for classifying the files
load('./FilesNeeded/Subjects.mat'); % A file (string array) with the subject names. It is supposed to be in a sub-folder of the current directory called FilesNeeded
load('./FilesNeeded/Times.mat');  % A file (string array) with the times [Variable1]. It is supposed to be in a sub-folder of the current directory called FilesNeeded
load('./FilesNeeded/Conditions.mat'); % A file (string array) with the conditions [Variable2]. It is supposed to be in a sub-folder of the current directory called FilesNeeded

%Define the saving folder
savingFolder = './MusicPlast_WithinBetween_MI/'; % change accordingly

%Select files on the basis of subject
for s=1:length(Subjects)   
subjectName = Subjects(s);
%Select files on the basis of condition 1
for t=1:length(Times) 
TimeName = Times(t);    
%Select files on the basis of condition 2
for c=1:length(Conditions)
conditionName = Conditions(c);
%Select all files satisfying condition1 condition2 subject selection 
selectedfiles = filedetails.Filename(filedetails.Subject == subjectName & filedetails.Condition == conditionName & filedetails.Time == TimeName);
%Select each of the files, one by one...
for f=1%:length(selectedfiles)
WorkFile = importdata(selectedfiles(f));
   
%Start MI estimation
WorkVariable = WorkFile.TF; % Name of the variable to work with

%Define parameters.
numSamples = 4001; %Number of samples in your trial
minDelay = 16; %Minimum delay for MI
maxDelay = 16; %Max delay for MI
numberOfFreqBands = 4; %How many bands are you interested in?
numSources = 360; %Number of sournces in your atlas

%Caclulate initial tables based on the parameters defined above
numSamplesMinusMinDelay = numSamples-minDelay; 
mi=zeros(1,numSamples);
miAtlasSingleSource=zeros(1,numSources);
miAtlasAllSources=zeros(numSources,numSources);
miAtlasAllPhaseBands=zeros((numSources*numberOfFreqBands),(numSources*numberOfFreqBands));
atlasLength=[1:numSources:(numSources*(numberOfFreqBands+1))];atlasLength(numberOfFreqBands+1)=(numSources*(numberOfFreqBands));

%Estimate Within and Between Bands Time delayed Mutual Info Of Phase
for pOne=1:numberOfFreqBands%for all phase bands
    for pTwo=1:numberOfFreqBands
        for k=1:numSources%for all sources as 1st part of comparison
            for i=1:numSources%for all sources as 2nd part of comparison
                %search mutual information of phase
                for d=minDelay:maxDelay%minDelay:numSamplesMinusMinDelay %search across all delays (from minDelay+)// or Choose 16ms based on https://doi.org/10.1016/j.neuroimage.2013.11.047
                    mi(1,d)=mi_gg((copnorm(((circshift((WorkVariable(k,:,pOne)),d,2))'))), (copnorm(((WorkVariable(i,:,pTwo)))')));
                  %  mi(1,1:minDelay)=0; mi(1,numSamplesMinusMinDelay:numSamples)=0; %zero out 1st and last
                    miMax=max(mi);%take the max of all delays
                end
                miAtlasSingleSource(1,i)=miMax;
            end
            miAtlasAllSources(k,:)=miAtlasSingleSource;
            miAtlasAllSourcesNorm = normalize(miAtlasAllSources);
        end
        miAtlasAllPhaseBands(atlasLength(pOne):((atlasLength(pOne)+(numSources-1))),atlasLength(pTwo):((atlasLength(pTwo)+(numSources-1))))=miAtlasAllSourcesNorm;
    end
end

%Create the outcome table for this file (trial)
trialHypernetwork = miAtlasAllPhaseBands;

% End MI estimation - Add each trial to a subject matrix.
SubjectTimeConditionTrialMatrix=zeros((numSources*numberOfFreqBands),(numSources*numberOfFreqBands),f);
SubjectTimeConditionTrialMatrix(:,:,f)= trialHypernetwork;

end
% Average MI files of each trial to create one hypernetwork
SubjectTimeConditionMatrix= mean(SubjectTimeConditionTrialMatrix,3);


%Save output file for subect-condition1-condition2
combinedStr = strcat(subjectName,'_',TimeName,'_',conditionName); %replace name!!!
combinedStr2 = strcat(savingFolder,combinedStr);
save(combinedStr2,'SubjectTimeConditionMatrix');   
end
end
end


% Function to be used for finding the correct files. It creates a table
% that has the characteristics of the file names [subject names, condition1, condition2, complete file name]
function filedetails = getdetails(filelist) %filelist: a structure array generated in the begining of the script
    filenames = vertcat(filelist.name); %convert the list of filename into a char array
    %temporary variables to create the table. Using string to make comparison easier
    Subject = string(filenames(:, 1:5)); % Subject names
    Time = string(filenames(:, 7:10)); % Post or Prex
    Condition = string(filenames(:, 12:17)); % Standx / AudDev / AvDevx / VisDev
    Filename = string({filelist.name}'); %complete file name
    %create the table
    filedetails = table(Subject, Time, Condition, Filename);
end