files = [{'DT3_rLS_rCA1_20150930_804um_576um'},...
    {'DT3_rLS_rCA1_20151002_948um_648um'},...
    {'DT3_rLS_rCA1_20151005_1380um_864um'},...
    {'DT3_rLS_rCA1_20151006_1524um_1224um'},...
    {'DT3_rLS_rCA1_20151006_1668um_1296um'},...
    {'DT3_rLS_rCA1_20151007_1740um_1296um'},...
    {'DT3_rLS_rCA1_20151008_1812um_1368um'},...
    {'DT3_rLS_rCA1_20151008_1956um_1440um'},...
    {'DT3_rLS_rCA1_20151009_2100um_1512um'},...
    ];

for f = length(files):-1:1;
    cd(files{f})
    xml = dir('*xml');
    SetCurrentSession(xml.name)
    units = GetUnits();
    units(units(:,1)<9,:)=[];
    for n = 1:size(units,1);
        spikes{n} = GetSpikeTimes(units(n,:));
    end
%     amp_corr =[];
%     amp_corr_flip = [];
%     if exist([files(f).name '-states.mat'])
%         load([files(f).name '-states.mat'])
%     end
%      
%     lfp = GetLFP([13 50]);
%     lfp = GetLFP([2 5 12 15 18 23 28 29]);
   lfp=GetLFP([29 60]);
%    if exist('ripples.rLS.evt')
%        delete('ripples.rLS.evt')
%    end
%    [b a] = butter(3,[140/625 190/625],'bandpass');
%    [ripples,stdev,noise] = FindRipples([lfp(:,1) filtfilt(b,a,lfp(:,2))],...
%        'durations',[25 120]); %,...
% %        'noise',[lfp(:,1) filtfilt(b,a,lfp(:,3))]);
%    SaveRippleEvents('ripples.rLS.evt',ripples,29);
%     save('ripples.mat','ripples','stdev','noise')
%     clear ripples stdev noise

%% waveforms
count=1;
for t = 1:units(end,1)
    u = find(units(:,1)==t);
    for i=1:length(u)
        temp = GetSpikeWaveforms(units(u(i),:));
        waveforms{t}(i,:,:)=squeeze(mean(temp));
        
        getfeatures = LoadSpikeFeatures([files{f} '.clu.' num2str(t)],20000);
        [L(count), Lratio(count), df(count)] = L_Ratio(getfeatures(:,4:end), find(getfeatures(:,3)==units(u(i),2)));
        IsolDist(count) = IsolationDistance(getfeatures(:,4:end), find(getfeatures(:,3)==units(u(i),2)));
        count=1+count;
    end
end
save('waveforms.mat','waveforms','L','Lratio','df','IsolDist','units','-v7.3')
clear ii count clusters kwx f u shank waveforms L Lratio IsolDist df

%% phaselocking
    if ~exist(['phaselocking.mat'])
        lfp = GetLFP([2 5 13 18 14 19 25 27 36 40 41 50 55 60 59]);
        n=1;
        for ch = 1:size(lfp,2)-1
            cells = find(units(:,1)==ch);
            if ~isempty(cells)
                
            for h = 1:100
               [b a] = butter(3,[h/625 (h*2+4)/625],'bandpass'); 
               lfp_filt = (filtfilt(b,a,lfp(:,ch+1)));
               phase = angle((hilbert(lfp_filt)));
               thresh = mean(fastrms(lfp_filt,round(1000/h))) + std(fastrms(lfp_filt,round(1000/h)));
               phase(fastrms(lfp_filt,round(1000/h))<thresh)=nan;
               
               for c = 1:length(cells)                  
                       histo{n,h} = phase(round(spikes{n}*1250));
                       histo{n,h}(isnan(histo{n,h})) = [];
                       if length(histo{n,h}) > 50
                       [phl_p(n,h) phl_z(n,h)] = circ_rtest(histo{n,h});
                       end
                       n=1+n;
                end
              
               n=1;
            end 
             end 
            ch
        end
        amp_corr = mean(amp_corr,3);

        save('phaselocking_power_thresh.mat','histo','phl*','-v7.3')

        clear amp* phl* histo
    end
%     

%% all assemblies
% if ~exist(['assemblies.mat'])
%     pair=[]
%     for i=1:length(units)
%         s(i) = length(spikes{i});
%     end
%     for i=1:length(units)
%         for j=1:length(units)
%     %         if units(i,1) < 9 && units(j,1) > 8  % implement for cross
%     %         regions only
%     %             if s(i) > median(s)-median(s)/2  && s(j) > median(s)  % sets a threshold on firing rate
%                     if i ~= j 
%                     pair = [pair;i j];
%                     end
%     %             end
%     %         end
%         end
%     end
% 
%     sp = zeros(length(units),round(length(lfp)/1250*1000));
%     for i=1:length(units)
%         sp(i,round(spikes{i}*1000))=1;
%     end
%     c=1;
%     for i=1:50000:length(sp)-50000
%         spp(:,c,:) = single(sp(:,i:i+50000));
%         c=1+c;
%     end
%     clear sp spikes lfp
%     u = unique(pair(:,1));
%     c=1;
%     for i=1:length(u)
%         f = find(pair(:,1)==u(i));
%         tic
%         disp(['running assembly analysis on ' xml.name ' on ' ...
%             num2str(length(f)) ' cell pairs, of ' num2str(length(pair)) ...
%             ' total pairs'])
%         l = size(spp,2);
%         [d dc] = cell_assembly(spp(:,1:round(l/2),:),0:150,pair(f,:));
%         [dd ddc] = cell_assembly(spp(:,round(l/2)+1:end,:),0:150,pair(f,:));
% 
%         for j =1:length(f)
%             dev(:,c,:) = squeeze(mean([mean(d(:,j,:),3) mean(dd(:,j,:),3)],2));
%             devControl(:,c,:) = squeeze(mean([mean(dc(:,j,:),3) mean(ddc(:,j,:),3)],2));
%             c=1+c;
%         end
%         save('assemblies.mat','dev*','pair','units')
%         toc
%     end
%     clear spp
%     save('assemblies.mat')
%     clear dev* pair 
% end
%     %% head direction
% % % 
% if ~exist('hd.mat')
% b = dir('Basle*.led');
% pos = importdata(b.name);
% frames = GetLFP(66);
% [pos] = alignPos(pos,frames);
% 
% for p=1:length(pos)
%     angl(p) = atan2(pos(p,2)-pos(p,4),pos(p,1)-pos(p,3));
%     if any(pos(p,1:4)==-1)
%         angl(p)=nan;
%     end
% end
% for i=1:length(units)
%     hd{i} = [];
% end
% 
% for i=1:length(units)
%     for j=1:length(spikes{i})
%     [a b] = min(abs(spikes{i}(j)-pos(:,5)));
%     if a < .034
%         hd{i} = [hd{i}; angl(b)];
%     end
%     end
% end
% 
% for i=1:length(units)
%     head_direction(i,:) = (hist([-pi; hd{i}; pi],36)./hist(angl,36));
% end
% save('hd.mat','-v7.3','head_direction','hd','angl','units')
% clear head_direction hd angl 
% end
%  %% home cage vs track 
% % % 
% if ~exist(['state_rate.mat'])
% b = dir('Basle*.led');
% pos = importdata(b.name);
% frames = GetLFP(66);
% [pos] = alignPos(pos,frames);
% 
% for i =1:length(units)
%    f = find(spikes{i}<pos(1,5));
%    ff = find(spikes{i}>pos(end,5));
%    
%    rate_homecage(i) = (length(f)+length(ff))./(length(find(frames(:,1)<pos(1,5)))+...
%        length(find(frames(:,1)>pos(end,5))))*1250;
%    rate_track(i) = (length(spikes{i})-(length(f)+length(ff)))./...
%        (length(intersect(find(frames(:,1)>pos(1,5)),...
%        find(frames(:,1)<pos(end,5)))))*1250;
% end
% 
% save('state_rate.mat','-v7.3','rate*')
% clear rate*
% end

%% power/rate correlations
%     if ~exist(['amp_rate_corrs.mat'])
%         for ch = 1:size(lfp,2)-1
%             amp_corr_temp = [];
%             amp_corr_flip_temp = [];
%             cells = find(units(:,1)==ch);         
%                   if ~isempty(cells)
%                       prate = zeros(length(cells),length(lfp));
%                       for c =1:length(cells)
%                         prate(c,round(spikes{cells(c)}*1250))=1;
%                         prate(c,:) = fastrms(prate(c,:),40);
%                       end
%                    f=1;
%                    cr = 1:1000000:length(lfp)-1000001;
%                    cr = cr(randperm(length(cr)));
%                    for chunk =  1:1000000:length(lfp)-1000001
%                       [ freqs, tt, spec ] = WaveSpec(lfp(chunk:chunk+1000000,ch+1),[1 200],100,5,1/1250,'log');
%                       spec = single(abs(spec));
%                       [ freqs, tt, specControl ] = WaveSpec(lfp(cr(f):cr(f)+1000000,ch+1),[1 200],100,5,1/1250,'log');
%                       specControl = single(abs(specControl));
%                     for c = 1:length(cells)
%                       for h = 1:100
%                           temp = corrcoef(abs(spec(h,200:1000000-200)),prate(c,chunk+200:chunk+1000000-200));
%                           amp_corr_temp(c,h,f) = temp(2);
%                           
%                           temp = corrcoef(abs(specControl(h,200:1000000-200)),prate(c,chunk+200:chunk+1000000-200));
%                           amp_corr_flip_temp(c,h,f) = temp(2);
%                           
%                       end
%                     end
%                     f=1+f;
%                    end
%                    end 
%             amp_corr = [amp_corr; amp_corr_temp];
%             amp_corr_flip = [amp_corr_flip; amp_corr_flip_temp];
%             ch
%             
%         end
%         amp_corr = nanmean(amp_corr,3);
%         amp_corr_flip = nanmean(amp_corr_flip,3);
% %         amp_corr_flip_std = nanstd(amp_corr_flip_temp,3);
% 
%         save('amp_rate_corrs.mat','amp_corr*','freqs','-v7.3')
% 
%         clear amp* 
%     end

    %% GTE
    
% %     
% if ~exist(['gte.mat'])
%     t = length(spikes);
%     sp = spikes;
%     sp{t+1} = [];
%     sp{t+2} = [t length(lfp)./1250];
%     [gte ci, all_te] = ASDFTE(sp,1:40);
%     clear sp t;
%     save('gte.mat','gte','ci','all_te','-v7.3')
% end

%% ripple modulation

% if ~exist(['ripple_modulation.mat'])
%     load('ripples.mat')
%     ripple_mod = zeros(size(units,1),1251);
%     [b a] = butter(3,[140/625 190/625],'bandpass'); 
%     n = 1;
%     for c = 1:length(units)
%         for i = 1:length(ripples)
% %             ll(n,:) = ll(n,:) + lfp_filt(locs_phase(i)-625:locs_phase(i)+625)';
% 
%         ff = find(spikes{n}>ripples(i,2)-.625);
%         fff = find(spikes{n}<ripples(i,2)+.625);
%         
%         ripple_mod(n,round((spikes{n}(intersect(ff,fff))-ripples(i,2))*1000)+626) = ...
%             ripple_mod(n,round((spikes{n}(intersect(ff,fff))-ripples(i,2))*1000)+626)+1;
%         end
%         numRipples(n) = length(ripples);
%         n=1+n;
%     end
% save('ripple_modulation.mat','ripple_mod','units','numRipples','-v7.3')
% clear avg* ll ripple* pks* locs* lfp_filt spec thresh numRipples
% end

%%     LFP phase and amplitude analysis
%  if ~exist(['lfp_phase_amplitude.mat'])
%      [b a] = butter(3,[4/625 9/625],'bandpass'); % filter bounds for theta
%      phase_amplitude = zeros(size(lfp,2)-1, 100,2501);
%      amplitude_amplitude = zeros(size(lfp,2)-1, 100,100);
%      phase_phase = zeros(size(lfp,2)-1, 99,361,361);
%      theta_wave = zeros(size(lfp,2)-1,2501);
%      
%      for ch = 1:size(lfp,2)-1       
%          for chunk = 1:1000000:length(lfp) - 999999
%              [ freqs, tt, spec ] = WaveSpec(lfp(chunk:chunk+1000000,ch+1),[1 200],100,5,1/1250,'lin');    
%              
%              % phase/amp
%              lfp_filt = (filtfilt(b,a,lfp(chunk:chunk+1000000,ch+1)));
%              [pks locs] = findpeaks(lfp_filt);
%              while locs(1) < 1251
%                  pks(1) = [];
%                  locs(1) = [];
%              end
%              while locs(end) + 1251 > length(spec)
%                  pks(end) = [];
%                  locs(end) = [];
%              end
%              for i=1:length(pks)
%                  phase_amplitude(ch,:,:) = squeeze(phase_amplitude(ch,:,:)) + abs(spec(:,locs(i)-1250:locs(i)+1250));
%                  theta_wave(ch,:) = theta_wave(ch,:) + lfp_filt(locs(i)-1250:locs(i)+1250)';
%              end
%  
%              % amp/amp
%              for f1 = 1:100
%                  for f2 = 1:100
%                      temp = corrcoef(abs(spec(f1,:)),abs(spec(f2,:)));    
%                      amplitude_amplitude(ch,f1,f2) = amplitude_amplitude(ch,f1,f2) + temp(2);
%                  end
%              end
% % 
%              % phase/phase
%              for f1 = 1:99
%                   [bb aa] = butter(3,[freqs(f1)/625 (freqs(f1+1)+1)/625],'bandpass');
%                   filt = (filtfilt(bb,aa,lfp(chunk:chunk+1000000,ch+1)));
%                   phase_theta = angle(hilbert(lfp_filt));
%                   phase_other = angle(hilbert(filt));
%                   for i = 1:length(phase_theta)
%                      t = round(rad2deg(phase_theta(i))+181);
%                      o = round(rad2deg(phase_other(i))+181);
%                      phase_phase(ch,f1,t,o) = phase_phase(ch,f1,t,o) + 1; 
%                   end
%                   ranges(f1,:) = [freqs(f1) freqs(f1+1)];
%              end
%          end            
%      end        
%   
%  save('lfp_phase_amplitude.mat','*phase*','ranges','*amplitude*','-v7.3')
%  end  
%  clear *phase* *amplitude* o t f1 f2 temp b a bb aa filt lfp_filt pks locs

 
 %% cross region analysis here
 
%  
%%    event trees 
% if ~exist(['eventtree.mat'])
%     [tree counts] = EventTree(spikes,.03,5,4);
% %     
% %     
%     save('eventtree.mat','tree','counts','-v7.3')
% %     
% end
    
    
    %% ACG analysis
%     if ~exist('acg.mat')
% frames = GetLFP(66);
% b = dir('Basle*.led');
% pos = importdata(b.name);
% [pos] = alignPos(pos,frames);
% 
%         sp = single(zeros(size(units,1),round(length(frames)./1250*1000)));
%         for i =1:size(units,1)
%            sp(i,round(spikes{i}*1000)) =1;
%         end
%         for i=1:size(units,1)
%             acg(i,:) = crosscorr(sp(i,:),sp(i,:),1000);
%             disp([num2str(((i-1)*size(units,1))./size(units,1)) ' percent done...'])
%         end
%         sp(:,1:round(pos(1,5)*1000))=0;
%         sp(:,round(pos(end,5)*1000):end)=0;
%         for i=1:size(units,1)
%             acg_behav(i,:) = crosscorr(sp(i,:),sp(i,:),1000);
%             disp([num2str(((i-1)*size(units,1))./size(units,1)) ' percent done...'])
%         end
%         
%         save('acg.mat','acg*','-v7.3')
%         clear sp acg*
%     end
% %     
%     %% CCG analysis
%     if ~exist('ccg.mat')
%         sp = single(zeros(size(units,1),round(length(lfp)./1250*1000)));
%         for i =1:size(units,1)
%            sp(i,round(spikes{i}*1000)) =1;
%         end
%         for i=1:size(units,1)
%             for j = i:size(units,1)
%                 ccg{i,j} = [];
%                 if units(i,1) ~= units(j,1)
% %                   if i==j
%                     ccg{i,j} = crosscorr(sp(i,:),sp(j,:),250);
%                 end
%             end
%             disp([num2str(((i-1)*size(units,1))./size(units,1)) ' percent done...'])
%         end
%         save('ccg.mat','ccg','-v7.3')
%         clear sp ccg
%     end
% %     
    cd ..
end
