SetCurrentSession
dirName = pwd;
units = GetUnits;
%units = units(units(:,1) ==2,:);
spikes  =GetSpikeTimes(units,'output','full');
spikestimes = spikes(:,1);
[uniqueCells,~,unitID] = unique(spikes(:,2:3),'rows');
spikesIDs = [spikes(:,2:3) unitID];
uniqueCells = sortrows(uniqueCells);
[spikestimes,b] = sort(spikestimes);
spikesIDs = spikesIDs(b,:);



fils  = getAllExtFiles(dirName,'ait',1);
stims = LoadEvents(fils{1});
sh = nan(length(stims.time/2),1);
for j = 1:2
   
    sti{j} = unique([stims.time(cellfun(@(a) any(regexp(a,'on')) & ...
        any(regexp(a,num2str(j))),stims.description)) ...
        stims.time(cellfun(@(a) any(regexp(a,'off')) & ...
        any(regexp(a,num2str(j))),stims.description))],'rows');
    
end

%categorize stimulations

fils  = getAllExtFiles(dirName,'aip',1);
stims1 = LoadEvents(fils{1});
temp = regexp(stims1.description,'\w*\d*','match');
pulse = cellfun(@(a) str2num(a{3}),temp);
sh = cellfun(@(a) str2num(a{2}),temp);

%pul1 = red light = PV stim
%pul2 = blue light = camkii stim

pul1 = pulse(sh==1);
pul2 = pulse(sh==2);

[a,b,c]  = histcounts(pul1,10);

kp1 = ismember(c,find(a>1000));
st1 = sti{1}(kp1,:);
cat1 = c(kp1);
[a,b,cat1] = unique(cat1,'sorted');
[a,b,c]  = histcounts(pul2,10);

kp2 = ismember(c,find(a>1000));
st2 = sti{2}(kp2,:);
cat2 = c(kp2);
[a,b,cat2] = unique(cat2,'sorted');

%find overlapping
[status,interval] = InIntervals(st1(:,1),st2);


%deal with overlaps

dual1_ts = st1(status);
dual1_cat = cat1(status);
dual2_ts = st2(interval(interval~=0));
dual2_cat = cat2(interval(interval~=0));
[a,b,dual_cat] = unique([dual1_cat dual2_cat],'rows');

single1_ts = st1(~status);
single1_cat = cat1(~status);


%%

popmat1 = nan(501,size(units,1),9);
%%
for i = 1:max(dual_cat)
    
   ts =  dual2_ts(dual_cat==i);
   eps = [ts-.1 ts+.5];
   for c=  1:size(units,1)
      
       spks = spikestimes(spikesIDs(:,3) ==c);
       spks = cellfun(@(a) linearize(histc(spks(spks>a(1) & spks<a(2))-a(1)-.1,-.1:.001:.4))',num2cell(eps,2), ...
           'uni',0);
       popmat1(:,c,i) = (nanmean(cell2mat(spks))) * 1000;
   end
    
end
%%
for c = [12]%18:size(popmat1,2)
figure
ax = tight_subplot(3,3);

for i = 1:9
    axes(ax(i))
    
     hold on
    plot([0 0],[0 60],'r')
     plot([0.3 0.3],[0 60],'r')
    
      plot([0.1 0.1],[0 60],'b')
     plot([0.15 0.15],[0 60],'b')
     
    bar( -.1:.001:.4,popmat1(:,c,i) )
    mx =max(max(popmat1(:,c,:),[],3));
    if mx>0
   ylim([0 mx+.1*mx])
    end
    xlim([-.1 .4])
end
mtit(num2str(uniqueCells(c,:)))
uiwait
end


%%



popmat = nan(501,3,size(units,1));
%popmat = nan(1001,3,size(units,1));
for i =  1:max(single1_cat)
    
   ts =  single1_ts(single1_cat==i);
   eps = [ts-.1 ts+.5];
   for c=  [6 11]%1:size(units,1)
      
       spks = spikestimes(spikesIDs(:,3) ==c);
       spks = cellfun(@(a) linearize(histc(spks(spks>a(1) & spks<a(2))-a(1),-.1:.001:.4))',num2cell(eps,2), ...
           'uni',0);
       
       
       popmat(:,i,c) = (nanmean(cell2mat(spks))) * 1000;
       %spks = cell2mat(spks);
     %  [p(c,i),h(c,i)] = signtest(mean(spks(:,1:100),2),mean(spks(:,201:300),2));
   end
    
end

%%
for c = 1:70
    
figure
ax = tight_subplot(1,3);
for i = 1:3
    axes(ax(i))
    
     hold on
    %plot([0 0],[0 60],'r')
     %plot([0.3 0.3],[0 60],'r')
    
      plot([0 0],[0 80],'b')
    
     
    bar( -.1:.0005:.4,popmat(:,i,c) )
    mx = max(max(popmat(:,:,c)));
    mx = mx+.1*mx;
     plot([0.05 0.05],[0 mx],'b')
  ylim([0 mx])
    xlim([-.1 .1])
end
uiwait
end
%%
bins = [-1:.5:4];
base = squeeze(nanmean(popmat(1:200,:,:)));
PYR = nanmean(base,1)<3;
evoked = squeeze(nanmean(popmat(201:300,:,:)));
gain = (evoked-base)./base;
distance = abs(units(:,1)-2);
col = linspecer(3,'jet');
ax = tight_subplot(1,3);

for j  = 1:3
    
    axes(ax(j))
mat = [];
for i = 0:2
    
    
   mat(i+1,:) = histc(gain(j,distance==i & PYR'),bins)/sum(histc(gain(j,distance==i & PYR'),bins));
end

 bar(bins,100*mat')
 if j==1
    legend('0um','200um','400um')
     ylabel('Percentage neurons')
 end
    xlabel('Percentage gain')
   
end
%%

global DATA
XMLfile = ['/media/sam/SM4/SeqGen/DS6/20161201/20161201.xml'];
[xml, rxml] = LoadXml(XMLfile);
ch_rip = [23 18 1 5];

[spks,b]  = sort(DATA.spikes(:,1));
[status,interval] = InIntervals(spks,sortrows(MergeEpochs2(cell2mat(sti'))));

st = status;
st(b) = status;
status = st;

waveform_ev = MeanSpikeShape_epoch(XMLfile,DATA.spikes(:,2:3),ch_rip,status);
waveform_sp = MeanSpikeShape_epoch(XMLfile,DATA.spikes(:,2:3),ch_rip,~status);
ch = cell2mat({xml.AnatGrps.Channels}')';
wav_ch = cell2mat({waveform_sp.ch});

%%
wav_ch = cell2mat({waveform_sp.ch});
b=[];
for i = 1:length(wav_ch)
    
    [a(i),b(i)] = find(ch == wav_ch(i));
    
end

ax = tight_subplot(1,3);

[~,a] = histc(a,0:3:9)
for i = 1:3
    axes(ax(i))
    ok  =accumarray([a' b'],gain(i,:)',[],@nanmean,nan);
    ok(isnan(ok)) = -10;
imagesc([1:4],[0:60:160],ok,[-1 5])
set(gca,'xtick',0:1:4,'ytick',0:60:160)
colormap('jet')
colordata = colormap;
colordata(1,:) = [1 1 1];
colormap(colordata);

    colorbar
ylabel('Vertical distance')
xlabel('Horizontal distnace')
end

%%
kp = ismember(units,[2 15],'rows');
    figure('visible','off')
    ax = tight_subplot(1,2);
    axes(ax(1))
    for j = 1:8
        
       plot(waveform_sp(kp).u(j,:)+j*700,'k')
       hold on
    end
    axis off
    title('Spont')
     axes(ax(2))
    for j = 1:8
        
       plot(waveform_ev(kp).u(j,:)+j*700,'k')
       hold on
    end
    
  axis off
    
    
    title('Evoked')
    
    %%
    
    
    mono_res1 = MonoSynConvClick3 (spikesIDs,spikestimes,'plot',false,'epoch',cell2mat(sti'),'binSize',.001);
mono_res2 = MonoSynConvClick3 (spikesIDs,spikestimes,'plot',false,'epoch',excludeEpochs([0 max(spikestimes)],cell2mat(sti')),'binSize',.001);
    %%
    figure
    kp = ismember(units,[2 15],'rows');
ax = tight_subplot(1,2);
axes(ax(1))
bar(-.1:.001:.1,(mono_res2.ccgR(:,kp,kp)/mono_res2.n(kp))/.001)
xlim([-.02 .02])
  ylim([0 80])
title('Spont')
axes(ax(2))
bar(-.1:.001:.1,(mono_res1.ccgR(:,kp,kp)/mono_res1.n(kp))/.001)
xlim([-.02 .02])

ylim([0 80])
    
    title('Evoked')
    
    %%
    for i = 1:70
        
     r(i)  = corr(waveform_sp(i).u(:),waveform_ev(i).u(:));
        
    end
    
    
    %%
    bar(.5:.01:1,histc(r,.5:.01:1))
    
    xlabel('Waveform Correlation')
    ylabel('Count')
    xlim([.6 1])
    
    
    %%
    
    %co stim
  stim = sortrows(cell2mat(sti'));
  stim = [stim(:,1)-.1 stim(:,2)+.1];
  [status, interval] = InIntervals(spikestimes,stim);
  
  
spk1 = spikestimes(spikesIDs(:,3)==20 & ~status);
spk2 = spikestimes(spikesIDs(:,3)==23& ~status);

times = [spk1;spk2];
groups = [ones(length(spk1),1);2*ones(length(spk2),1)];
[times,b] = sort(times);
groups = groups(b);


%%
close all
[ccg2,n1,t1] = CCGTime(times,groups,'binSize',.0005,'duration',.2,'win',1000,'win_inc',1000);
ok = squeeze(ccg2(:,1,2,:));
nn = squeeze(n1(:,1,2,:));
tr =[];
pr = [];
for i = 1:size(ok,2)
    [tr(i),pr(i,:),prob_uncor,pred] = GetTransProb(ok(:,i),nn(1,i),.0005);
end


pp = cumsum(histc(dual2_ts,t1{2}));

plot(log10(pp),tr,'x')
%%
  ts =  dual2_ts(dual_cat==9);
  
  for i = 1:size(ts)
      
      
      pre = [ts(i)-1000 ts(i)];
       post = [ts(i) ts(i)+1000];
       
       s1 = spk1(spk1>pre(1) & spk1<pre(2)); 
       s2 = spk2(spk2>pre(1) & spk2<pre(2));
       rawCCG = CrossCorr(s1,s2,.0005,501);
       [trans_pre(i),prob_pre(i,:),prob_uncor,pred] = GetTransProb(rawCCG,length(s1),.0005);
       
        s1 = spk1(spk1>post(1) & spk1<post(2)); 
       s2 = spk2(spk2>post(1) & spk2<post(2));
       rawCCG = CrossCorr(s1,s2,.0005,501);
       [trans_post(i),prob_post(i,:),prob_uncor,pred] = GetTransProb(rawCCG,length(s1),.0005);
       
       
  end
  
  
  %%
 
  %%
  
  
    %not co stim
    tr =[];prob_uncor=[];
pr = [];
  gd_pair =   pre_post(pre_post(:,2)==12,:);
  gd_pair = gd_pair(1:end-4,:);
  spk2 = spikestimes(spikesIDs(:,3)==12& ~status);
  
    for cc = 1:size(gd_pair,1)
  stim = sortrows(cell2mat(sti'));
  stim = [stim(:,1)-.1 stim(:,2)+.1];
  [status, interval] = InIntervals(spikestimes,stim);
  
  
spk1 = spikestimes(spikesIDs(:,3)==gd_pair(cc,1) & ~status);


times = [spk1;spk2];
groups = [ones(length(spk1),1);2*ones(length(spk2),1)];
[times,b] = sort(times);
groups = groups(b);


[ccg2,n1,t1] = CCGTime(times,groups,'binSize',.001,'duration',.2,'win',1000,'win_inc',1000);
ok = squeeze(ccg2(:,1,2,:));
nn = squeeze(n1(:,1,2,:));

for i = 1:size(ok,2)
    [tr(i,cc),pr(i,:,cc),prob_uncor(i,:,cc),pred] = GetTransProb(ok(:,i),nn(1,i),.001);
end
    end

pp = cumsum(histc(dual2_ts,t1{2}));

plot((pp),tr)

%%

 close all
  y = nanmedian((tr-repmat(tr(1,:),size(tr,1),1)),2);
   z = SEM((tr-repmat(tr(1,:),size(tr,1),1)),2);
   plotShadedError(pp(1:end-1),y(1:end-1),z(1:end-1),[0 0 0],.25)
   hold on
  plot(pp(1:end-1),y(1:end-1),'k','linewidth',3)
  figure
  bar(t1{1},pr(3,:,11))
  ylim([-.01 .06])
    xlim([-.1 .1])
    figure
  bar(t1{1},pr(22,:,11))
    ylim([-.01 .06])
     xlim([-.1 .1])
%%
  ts =  dual2_ts(dual_cat==9);
  
  for i = 1:size(ts)
      
      
      pre = [ts(i)-1000 ts(i)];
       post = [ts(i) ts(i)+1000];
       
       s1 = spk1(spk1>pre(1) & spk1<pre(2)); 
       s2 = spk2(spk2>pre(1) & spk2<pre(2));
       rawCCG = CrossCorr(s1,s2,.0005,501);
       [trans_pre(i),prob_pre(i,:),prob_uncor,pred] = GetTransProb(rawCCG,length(s1),.0005);
       
        s1 = spk1(spk1>post(1) & spk1<post(2)); 
       s2 = spk2(spk2>post(1) & spk2<post(2));
       rawCCG = CrossCorr(s1,s2,.0005,501);
       [trans_post(i),prob_post(i,:),prob_uncor,pred] = GetTransProb(rawCCG,length(s1),.0005);
       
       
  end
  

  %%
  maxw = [];
  for i = 1:4
  filename = ['/media/sam/SM4/SeqGen/DS6/20161201/20161201.spk.' num2str(i)];
  waveforms = LoadSpikeWaveforms(filename,8,32);
    w =   double(min(waveforms(:,:,16),[],2));
     filename = ['/media/sam/SM4/SeqGen/DS6/20161201/20161201.clu.' num2str(i)];
clu = load(filename);

w = w(clu(2:end)>1);

maxw = [maxw;w];
  end