global DATA
import CMBHOME.*
if ~isstruct(DATA)
    [filename,path] = uigetfile('*.xml','Please select a parameter file for this session');
    SetCurrentSession([path filename])
    [xml, rxml] = LoadXml(filename);
    
end
path = DATA.session.path;

sl = regexp(path,'/');
pardir = path(1:sl(end));
%get analogin

fils = getAllDatFiles(pardir);

fils = fils(cellfun(@(a) any(regexp(a,'auxiliary')),fils));
cmd = [];
if length(fils)>1
    for i = 1:length(fils)
    cmd = [cmd fils{i}];
    end
end
%concat all dat files, downsample and save


%allocate memory
d = nan(4,1250*60*60*4);siz = 0;
for i = 1:length(fils)
    if i ==1
        %first session only two analogin channnels turned on
    temp = LoadBinary(fils{i},'nChannels',2,'channels',[1:2]);
    temp = [temp temp];
    else
         temp = LoadBinary(fils{i},'nChannels',4,'channels',[1:4]);
    end
    temp = temp(1:16:end,:);
 
    d(:,siz+1:siz+size(temp,1)) = temp';
       siz = size(temp,1)+siz;
       len(i) = size(temp,1);
end
outfile = [DATA.session.basename '_CMB.mat'];
analogin = d(:,~any(isnan(d)));
ts = 1/1250:1/1250:length(d)/1250;
b_ts = 1/39.125: 1/39.125:max(ts);
x = nan(length(b_ts),1);
y = nan(length(b_ts),1);
hd = nan(length(b_ts),1);
event = [];
root = CMBHOME.Session('name', outfile, 'b_ts', b_ts,...
    'b_x',x, 'b_y', y, 'b_headdir', hd,...
    'event', {}, 'raw_pos', 1, 'raw_headdir', 1, 'date_created', now, ...
    'epoch', ts([1 end]), 'fs_video', 39.0625, 'path_raw_data', path);

%%



units = GetUnits;
spikes = GetSpikeTimes(units,'output','numbered');


for ii = 1:size(units,1)
    
    spk_ts = spikes(spikes(:,2)==ii,1);
    spk_i = time2ind(root.b_ts,spk_ts);
    spk_ts = spk_ts(spk_i>0);
    spk_i = spk_i(spk_i>0);
    
    root.spike(units(ii,1),units(ii,2)) = CMBHOME.Spike('i',spk_i,'ts',spk_ts,'tet_name',num2str(units(ii,:)));
end

%%

root.path_lfp = repmat({[path '/' filename]},128,1);

%%


%get lfp


fils = getAllExtFiles(pardir,'eeg');


%concat all dat files, downsample and save


%allocate memory
lfp = nan(1,length(analogin));siz = 0;
for i = 1:length(fils)
    
    temp = LoadBinary(fils{i},'nChannels',32,'channels',24);
   
 
    lfp(:,siz+1:siz+size(temp,1)) = temp';
       siz = size(temp,1)+siz;
end



%%



root.path_lfp = repmat({[path '/' filename]},5,1);


    
    fs = 1250;
    channel_name=['sh 1 ch 5'];
    
    root.b_lfp(1) = CMBHOME.LFP(lfp, ts, fs, channel_name);
    root.b_lfp(1) = root.b_lfp(1).AppendTheta;
    root.b_lfp(1) = root.b_lfp(1).AppendThetaAmplitude;
    root.b_lfp(1) = root.b_lfp(1).AppendThetaPhase;
    
    if ~isfield( root.lfp_manager,'theta_freq')
        root.lfp_manager.theta_freq=0;
    end
    root.b_lfp(1) = root.b_lfp(1).AppendThetaFreq;
    
    
    root.b_lfp(1).channel_name = channel_name;
    
    for i = 1:4
            channel_name=['analogin' num2str(i)];
         root.b_lfp(1+i) = CMBHOME.LFP(analogin(i,:), ts, fs, channel_name);
        
    end
end

%%
for ii = 1:size(root.cells,1)
    
    spk_lfp = Utils.SpkInds(ts,root.spike(root.cells(ii,1),root.cells(ii,2)).ts);
    root.spike(root.cells(ii,1),root.cells(ii,2)).i_lfp = spk_lfp;
end

idx = 1;

thres = [8000 4300 4300 4300];
for i = 2:5
    
    
    [~,pks] = findpeaks(root.b_lfp(i).signal,'MinPeakDistance',1250,'MinPeakHeight',thres(i-1));
    
    for j = 1:length(pks)
   ev{idx,1} =  ts(pks(j));
  ev{idx,2} =   ['anin ' num2str(i-1)];
  idx = idx+1;
end
end


[~,b]  = sort(cell2mat(ev(:,1)));

root.event = [ev(:,2) ev(:,1)];
save([path '/' outfile],'root','-v7.3')

%%
ts = root.Label2Time('anin 1');
binnedPop=populationMatrix(root,.5,1.5,200,ts,k');
k = gaussian2Dfilter([20 1],1);
kp = mean(mean(binnedPop,2),3)<5 & max(mean(binnedPop,3),[],2)>5;
[xx,yy] = meshgrid(root.cells(kp,1));
up = ones(size(xx));
up = triu(up,1);
dt = 1:10:length(states);
idx = 1; cor_sam = []; cor_dif = [];
for i = 1:10:dt(end)
binnedPop=populationMatrix(root,0,100,10000,i);
mat = nanzscore(binnedPop(kp,:),[],2);
a = corr(mat');
cor_sam(idx) = nanmean(a(up & xx==yy));
cor_dif(idx) = nanmean(a(up & xx~=yy));
idx= idx+1;
end
%%

%dt1 = 1:length(states);
%cor_sam1 = interp1(1:10:15500,cor_sam,1:length(states));
%cor_dif1 = interp1(1:10:15500,cor_dif,1:length(states));
r = []; p = [];
for i = 1:length(cor_sam) - 100
[temp1,temp2] = corrcoef(cor_sam(i:i+100),cor_dif(i:i+100),'rows','pairwise');
r(i) = temp1(1,2);
p(i) = temp2(1,2);
end

%%


figure
ax = tight_subplot(3,1);

axes(ax(1))
plot(states,'k','linewidth',3)
set(gca,'ydir','reverse','xticklabel',[])
xlim([0 length(states)])
hold on
plot([4203.208 4203.208], [-1 6],'k')
plot([10272.7872 10272.7872], [-1 6],'k')
ylim([-1 6])

axes(ax(2))
%bar(dt,cor_sam-cor_dif)
bar(dt(1:length(r)),r)
xlim([0 length(states)])
hold on
%plot([4203.208 4203.208], [-1 1],'k')
%plot([10272.7872 10272.7872], [-1 1],'k')

set(gca,'xticklabel',[])
axes(ax(3))

[hAx,hLine1,hLine2] = plotyy(dt(1:length(cor_sam)),cor_sam,dt(1:length(cor_sam)),cor_dif);
hLine1.Color = 'k';
hLine1.LineWidth = 2;
hLine2.Color = 'r';
hLine2.LineWidth = 2;

set(hAx(2),'fontsize',16)
hold on

plot([4203.208 4203.208], [-.02 .05],'k')
plot([10272.7872 10272.7872], [-.02 .05],'k')