
% TET-ANALYSIS AND DATABASE FOUNDATION
% MV-EBA 2018  

% Superfolder with the experiments got by tetModel
%cd('/Volumes/Enhorhinal/Dropbox/TETMODEL/tetModel/DATA/'); % Manu Macbook
cd('D:\Dropbox\TETMODEL\tetModel\DATA'); % ManuPC
%cd('C:\Users\estibaliz\Dropbox\LABORATIORIO 2015-2019 organised\LAB\2018\FEBRERO 18\MODELO CIRCUITO NEURONAL\MODELO RipTet\tetModel\DATA'); % TetPC
allFolder = genpath(pwd); allFolder = strsplit(allFolder,';');% get all folders
if length(allFolder)==1 % for apple systems
    allFolder = genpath(pwd); allFolder = strsplit(allFolder,':');% get all folders
end

database = []; waveAll = []; freq = []; exp = {};
for ii = 1:size(allFolder,2) % go for all folders
    if contains(allFolder{ii},'tetSim') % if folder contain an experiment
        disp(' ');
        fprintf('%i of up to %i... \n',ii, size(allFolder,2));
        disp(allFolder{ii});
        cd(allFolder{ii}); % go to the experiment folder
        try error %load('TetData.mat'); % if already exist
            disp('Loaded TetData.mat...');
        catch
            disp('Loading data...');
            try mv = csvread('mv.cvs');
                expName = strsplit(pwd,{'\','/'}); expName = expName{end};
                feat = csvread('feat.cvs',1); % load experimental data and description, 
                log = csvread('log.cvs',1); % column 1, cell type, column 2, tet
                spk = csvread('spk.cvs');
                
                % 0 GET FEATURES
                d(1) = feat(15); % Code, 0 is control, 1 is all, 2 is lead, 3 is exc and 4 is inh
                d(2) = feat(14); % Prob tet, from 0-1
                if feat(7) == 0 && feat(8) == 0 && feat(9) == 0 % if no synapsis
                    d(1) = 5; % For no synapsis/ control networks
                end
                dur = feat(2); % Duration in seconds

                exc_cells = find(log(:,1)==0 & log(:,2)<1); % all exc cells, non-considering tets
                inh_cells = find(log(:,1)==2 & log(:,2)<1); % all inh cells, non-considering tets
                lead_cells = find(log(:,1)==1 & log(:,2)<1); % all leaders cells, non-considering tets
                all_cells = find(log(:,1)<3 & log(:,2)<1); % all leaders cells, non-considering tets

                % 1 FIRING RATE AND AVERAGE MP
                disp('Computing firing rate...');
                d(3) = length(find(ismember(spk(1,:),exc_cells)))/length(exc_cells)/dur; % exc cells
                d(4) = length(find(ismember(spk(1,:),inh_cells)))/length(inh_cells)/dur; % inh cells     
                d(5) = length(find(ismember(spk(1,:),lead_cells)))/length(lead_cells)/dur; % lead cells
                d(6) = length(find(ismember(spk(1,:),all_cells)))/length(all_cells)/dur; % all cells

                % Average mp
                win = find(mv(:,end) > 0.05); % discarting 50 first milliseconds 
                d(7) = mean(mean(mv(win,exc_cells),2));
                d(11) = std(mean(mv(win,exc_cells),2));

                d(8) = mean(mean(mv(win,inh_cells),2));
                d(12) = std(mean(mv(win,inh_cells),2));

                d(9) = mean(mean(mv(win,lead_cells),2));
                d(13) = std(mean(mv(win,inh_cells),2));

                d(10) = mean(mean(mv(win,all_cells),2));
                d(14) = std(mean(mv(win,all_cells),2));

                % 2 NETWORK CORRELATION
                disp('Computing synchrony...'); % correlation defined as the mean correlation beetween vm
                % synchrony exc cells
                C = corr(mv(win,exc_cells)); % get correlation for mv Exc cells
                d(15) = mean(C(:)); % mean synch
                d(19) = median(C(:)); % median synch
                d(23) = std(C(:)); % std synch
                % inh cells
                C = corr(mv(win,inh_cells)); % get correlation for mv Ing cells
                d(16) = mean(C(:)); % mean synch
                d(20) = median(C(:)); % median synch
                d(24) = std(C(:)); % std synch
                % lead cells
                C = corr(mv(win,lead_cells)); % get correlation for mv Ing cells
                d(17) = mean(C(:)); % mean synch
                d(21) = median(C(:)); % median synch
                d(25) = std(C(:)); % std synch
                % all cells
                C = corr(mv(win,all_cells)); % get correlation for mv Ing cells
                d(18) = mean(C(:)); % mean synch
                d(22) = median(C(:)); % median synch
                d(26) = std(C(:)); % std synch

                % 3 SPECTRAL FEATURES
                disp('Computing spectral features');
                lfp = -mean(mv(:,1:end-1),2); lfp = lfp - mean(lfp); % estimate lfp as all membrane current average
                fs = 1/mean(diff(mv(:,end))); % calculate sampling frequency
                wav = wavelet_format(lfp(win), fs, [1 100],1,1); % compute wavelet

                w = mean(abs(wav.cfs),2); % storing wavelet analysis
                freq = wav.frequencies;
                % mean power
                d(27) = mean(w(freq>=28 & freq<=65)); % take spectrum values, gamma
                d(28) = mean(w(freq>=16 & freq<=27)); % beta
                d(29) = mean(w(freq>=6 & freq<=16)); % theta
                d(30) = mean(w(freq>=1 & freq<=6)); % delta
                % max power
                d(31) = max(w(freq>=28 & freq<=65)); % take spectrum values, gamma
                d(32) = max(w(freq>=16 & freq<=27)); % beta
                d(33) = max(w(freq>=6 & freq<=16)); % theta
                d(34) = max(w(freq>=1 & freq<=6)); % delta
                % normalize spectrum
                wN = w/sum(w);
                d(35) = max(wN(freq>=28 & freq<=65)); % take spectrum values, gamma
                d(36) = max(wN(freq>=16 & freq<=27)); % beta
                d(37) = max(wN(freq>=6 & freq<=16)); % theta
                d(38) = max(wN(freq>=1 & freq<=6)); % delta

                save('TetData.mat','d','expName','w','freq');
            catch
                disp('No mv! Skip!');
                d = []; expName = []; w = [];
            end
        end
        database = [database; d];
        exp{size(exp,2)+1} = expName;
        waveAll = [waveAll; w'];
        clear d expName w
    end
end

% saving database, waveAll and freq
cd(allFolder{1});
save(strcat('TetDatabase_',date,'.mat'), 'database', 'freq', 'waveAll','exp');
