function [PDC,Signif,S,H,COH,DTF,dDTF,pCOH,ep_times,E,freq] = MainPDC_spikes_times

% Handles PDC analysis for spike train-type data only.
% 1) Performs PCA on sets of spike data (or loads PCs directly)
% 2) Divides data into epochs
% 3) Calculates optimal model order (optional)
% 4) Creates MVAR model (taken from biosig!)
% 5) Checks for model stability (stationarity)
% 6) Calculates PDC and other measures (biosig)
% 7) Calculates PDC-siginificance level (Schelter et al. 2005)
% 8) Runs PDC validation tests
% 9) Plots results.
% NOTE: the different spike trains-groups MUST have the same binning in
% order to model them with MVAR models. (???)


% Perform PCA on groups of spike trains....................................
disp('Do you want to preprocess spike trains or are you going to use pre-processed data from a file?')
disp('Process spikes -> 1');                                                % Choose whether to process spike trains...
disp('Load processed spikes from a file -> 2');                             % or load the PCs directly.
work_type = input('');
if work_type ~= 1 && work_type ~= 2                                         % If none of the 2...
    error('You have to choose one of the two options');                     % then stop!
end

if work_type == 1                                                           % If process spike trains...
    [pr_spikes,PC1,f_N,maxtime,mintime,FileName] = Proc_spikes;                     % Call Prc_spikes to do the analysis and get PC1.
    %FileName=['Spikes_',FileName];
    count_rate=2*f_N;                                                       % The sampling rate is set as double the Nyquist frequency.
    clear pr_spikes
    
elseif work_type == 2
    [FileName,PathName] = uigetfile('*.mat','Select the Mat-file with the pre-processed data'); % If load file
    if isequal(FileName,0)
        error('You did not specify a file')
    else
        disp(['Opening file: ', fullfile(PathName, FileName)])
    end
    h = fullfile(PathName, FileName);
    load(h);                                                                % Load the file (to LOAD xreiazeai gia na anoixei to archive!)
    count_rate = archive.sampling_rate;                                     % and the necessary variables from the archive
    exp_time = archive.experim_time;
    mintime = exp_time(1); maxtime = exp_time(2);
end


% Divide data into epochs..................................................
bins = repmat([mintime:1/count_rate:maxtime]',1,2);                         % Create a bins array to run the SetEpochsPCA 
[dataep,ep_times,E] = SetEpochsPCA(PC1,bins);                               % Divide into epochs.

m = length(PC1);                                                            % Count number of variables

for i = 1:E
    mm = mean(dataep{i});                                                   % Calculate the data mean (columnwise)
    mm = repmat(mm,size(dataep{i},1),1);
    dataep{i} = dataep{i} - mm;                                             % and subtract them from the data
end

% Allocate memory for variables............................................
Ar = cell(1,E);                                                             % Ar = MVAR coefficients
Af = cell(1,E);                                                             % Af = Fourier transform of Ar
Sigma = cell(1,E);                                                          % Sigma = error covariance matrix
S = cell(1,E);                                                              % S = cross spectral density matrix
H = cell(1,E);                                                              % H = transfer matrix
PDC = cell(1,E);                                                            % PDC
Signif = cell(1,E);                                                         % Signif = PDC-significance level
COH = cell(1,E);                                                            % COH = coherency (take abs for magnitude)
DTF = cell(1,E);                                                            % DTF = directed transfer function
dDTF = cell(1,E);                                                           % dDTF = direct directed transfer function
pCOH = cell(1,E);                                                           % pCOH = partial coherence
sp = cell(1,E);                                                             % data spectral density (from Neurospec)


% Calculate optimal model order............................................
chooseP = input('Do you want to see the optimal MAR model order according to the SBC and FPE criteria? Y/N: ', 's');
if chooseP == 'Y' || chooseP == 'y'                                         % If user chooses to search for optimal model order...
    pmin = 1;                                                               % it will check model orders from pmin...
    pmax = input('Input the maximum model order you want to check (default 100): '); % to pmax.
    if isempty(pmax)
        pmax=100;
    end
    pmax=round(pmax);
    for i = 1:E                                                             % And for each epoch...
        [pbic, pfpe, logdpbic, logdpfpe] = ModelOrder2(dataep{i},pmin,pmax);% It will calculate optimal model order according to BIC and FPE
        disp(['BIC: The optimal model order for epoch ', num2str(i), ' is: ', num2str(pbic),' (logdp = ', num2str(logdpbic),')'])
        disp(['FPE: The optimal model order for epoch ', num2str(i), ' is: ', num2str(pfpe),' (logdp = ', num2str(logdpfpe),')'])
        disp(' ')
    end
end
    P = input('Input the MAR model order (i.e. history [e.g. 100 (ms)]): ');% Specify the model order explicitly
    disp(' ');
    P = round(P);

% Create MVAR model........................................................
disp('Specify frequency range for spectral analysis [e.g. 50 for range 0-50Hz]');
disp('Or input row matrix with specific frequencies')
disp(['Note: The maximum possible frequency is half the sampling rate = ', num2str(count_rate/2)]); % Inform user on the Nyquist freq.
freq = input('Frequencies = ');                                             % Specify the frequency range for the spectral analysis
disp(' ')                                                                   
freq = sort(freq);

for i=1:E                                                                   % For each epoch...
    disp(['Performing VAR analysis on data of epoch ' num2str(i) '.'])
    %disp(' ')
    [AR,RC,sigma] = mvar(dataep{i}, P, 2);                                  % Call mvar to perform a MAR analysis

    
% Check model stability....................................................
    disp('Performing VAR stability test')                                   % Test the stability of the VAR model...
    z = findroots(full(AR),m,P);                                            % by locating the eigenvalues of matrix A (Lutkepohl)
    if all(z <= 1)
        disp('Stability criterion is satisfied')
    else
        disp('Stability criterion is violated!');                           % If the model is unstable (eigenbvalues > 1) stop!
        disp('End of process!')
        return
    end
    %disp(' ')
    
    
% Calculate PDC etc........................................................
    B = eye(m);                                                             % B is the identity matrix
    A = [eye(m),-AR];                                                       % Bring Ar in the appropriate form...
    C = sigma(:,m*P+1:m*(P+1));                                             % and keep only the last M columns of sigma (noise covariance)

    if isscalar(freq)                                                       % If the frequency is just a number...
        freq=0:freq;                                                        % then set it to be the vector 0:freq
    end
    disp(['Calculating the power spectra and PDC coefficients of epoch ' num2str(i) '.'])
    %disp(' ')
    [AF,s,h,pdc,coh,dtf,pcoh,ddtf] = mvfreqz2(B,A,C,freq,count_rate);       % Call mvfreqz2 (modified by me) to calculate the desired measures

    disp(['Computing the PDC significance levels of epoch ' num2str(i) '.'])
    %disp(' ')
    signif=Significance(dataep{i},AF,C,freq,P,count_rate);                  % Calculate the significance level of PDC
    
    Ar{i} = AR;                                                             % Store them in cell arrays...
    Af{i} = AF;                                                             % where each cell corresponds to...
    Sigma{i} = sigma;                                                       % a specific epoch.
    S{i} = s;                                                                
    H{i} = h;                                                               
    PDC{i} = pdc;                                                           
    Signif{i} = signif;
    COH{i} = coh;                                                           
    DTF{i}= dtf;
    dDTF{i}= ddtf;
    pCOH{i} = pcoh;
    
   
    disp('Plotting the model accuracy');
    sp{i} = Spctr_comp(dataep{i}, S{i}, freq, count_rate, 7);              % Calculate and plot the data power-spectra
    disp(' ');
    
    dataep{i}=[];
    clear AR RC sigma s h pdc coh dtf pcoh ddtf AF signif A B C
end                                                                         % Repeat for the next epoch


% Run validity tests.......................................................
disp('Running tests of PDC computational efficiency')                       % Run tests for computational validation of PDC
PDCtest1(PDC,E)                                                             % check if PDC >= 0 & PDC <= 1
PDCtest2(PDC,E,m)                                                           % check if Sum(PDCi) = 1

archive.file_name = FileName;
archive.epochs = ep_times;
archive.model_order = P;
archive.sampling_rate = count_rate;
archive.frequencies = freq';

uisave({'archive','Ar','Af','Sigma','S','H','PDC','Signif','COH','DTF','dDTF','pCOH','sp','ep_times','freq'},['PDCresults_Spikes_', FileName]);
                                                                            % save in user-specified file

% Plot resuts..............................................................
PlotPDC(PDC,COH,S,Signif,ep_times,freq,E,m,1);                              % Call plotPDC to plot results