%get position for camkii7 10/14/17

basepath  ='/media/sam/SM4/SeqGen/camkii4/20160911';
basename = '20160911';

digitalin_ch = 1;

 %pos_inf =  GetPosBaslerBlinkingLED(basepath,[],basename,digitalin_ch);
load('position_info.mat')

del = false(length(pos_inf.in_eps),1);
for i = 1:length(cellfun(@range,pos_inf.len_ep))-1
    if pos_inf.in_eps(i) == pos_inf.in_eps(i+1)
 del(i+1) = true;
         pos_inf.ts_ep{i} = [pos_inf.ts_ep{i};pos_inf.ts_ep{i+1}];
         pos_inf.len_ep{i} = [pos_inf.len_ep{i};pos_inf.len_ep{i+1}];
    end
end
pos_inf.len_ep(del) = [];
pos_inf.ts_ep(del) = [];
pos_inf.in_eps(del) = [];
pos_inf.out_eps(del) = [];
%merge trials



    fils  = getAllExtFiles(basepath,'ait',1);
stims = LoadEvents(fils{1});

for j = 1:5
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

fils  = getAllExtFiles(basepath,'aip',1);
stims = LoadEvents(fils{1});
ok = regexp(stims.description,'\w*','match');
amp  =cellfun(@(a) str2num(a{3}),ok);
sh  =cellfun(@(a) str2num(a{2}),ok);
%%
SetCurrentSession
units = GetUnits;
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,:);

%%
ineps = cell2mat(cellfun(@(a) [a(1) a(end) ] ,pos_inf.ts_ep(pos_inf.in_eps),'uni',0));
outeps = cell2mat(cellfun(@(a) [a(1) a(end) ] ,pos_inf.ts_ep(pos_inf.out_eps),'uni',0));



    
    
    
   

t= pos_inf.ts;
%kp = InIntervals(sti{1}(:,1),[min(t) max(t)]);
len_out_st = [];
len_out_fin = [];
interval_out =[];

for i = 1:4
   [status,interval] = InIntervals(t,sti{i}(:,:));
   len = pos_inf.lin_pos;
   n1 = histoc(interval(interval>0),1:size(sti{i}(:,1),1));
    len_st = mat2cell( len(status),n1);
    
    
    len_out_st{i} = cellfun(@min,len_st(cellfun(@any,len_st)));
    len_out_fin{i} = cellfun(@max,len_st(cellfun(@any,len_st)));
    
    
    
       [~,interval_out{i}] = InIntervals(sti{i}(:,1),sortrows([outeps]));
   interval_out{i} = interval_out{i}(interval_out{i}>0);
   
   
   
    [status,interval] = InIntervals(t,sti{i}(:,:));
    
    
   len = pos_inf.lin_pos;
   n1 = histoc(interval(interval>0),1:size(sti{i}(:,:),1));
    len_st = mat2cell( len(status),n1);
    
      len_in_st{i} = cellfun(@min,len_st(cellfun(@any,len_st)));
    len_in_fin{i} = cellfun(@max,len_st(cellfun(@any,len_st)));
    
    
   [~,interval_in{i}] = InIntervals(sti{i}(:,1),sortrows([ineps]));
   interval_in{i} = interval_in{i}(interval_in{i}>0);
   
   
   
end
   
   %%
   
   
   
   
   
    stim.description = 'stim on S1,S2,S3,S4 directory before/after linear track. then gaussian at fixed pos on track';
    stim.ts = cell(4,1);
stim.amp = cell(4,1);
stim.ep = cell(4,1);
ineps = cell2mat(cellfun(@(a) [a(1) a(end) ] ,pos_inf.ts_ep(pos_inf.in_eps),'uni',0));
outeps = cell2mat(cellfun(@(a) [a(1) a(end) ] ,pos_inf.ts_ep(pos_inf.out_eps),'uni',0));



    
    t= pos_inf.ts;
    
   
fils  = getAllExtFiles(basepath,'aip',1);
stims = LoadEvents(fils{1});
ok = regexp(stims.description,'\w*','match');
sh  =cellfun(@(a) str2num(a{2}),ok);
amp  =cellfun(@(a) str2num(a{3}),ok);


for ii = 1:4
amp2 = amp(sh==ii);
pre2 = sti{ii}(:,1)<5600;
post2 = ( sti{ii}(:,1)>8600);



temp2 = amp2(post2|pre2);
[a,b] = histcounts(temp2,5);
[a,b2] = histc(temp2,b);




stim.ts{ii} = sti{ii};
stim.amp{ii} = nan(size(sti{ii},1),1);
stim.amp{ii}(pre2|post2) = b2;
stim.ep{ii} = cell(size(sti{ii},1),1);
stim.ep{ii}(pre2) = repmat({'pre'},sum(pre2),1);
stim.ep{ii}(post2) = repmat({'post'},sum(post2),1);
stim.ep{ii}(~post2& ~pre2) = repmat({'during'},sum(~post2& ~pre2),1);


end


save('stim_info.mat','stim')

    %%
    rate2 = nan(5,2,size(uniqueCells,1));
    rate3 = nan(5,2,size(uniqueCells,1));
    
for i = 1:size(uniqueCells,1)
   
    spks = spikestimes(spikesIDs(:,3)==i);
    [status,interval] = InIntervals(spks,eps2);
    
   
    
    n1 = histoc(interval(interval>0),1:size(eps2,1));
    
    p1 = mat2cell(spks(status),n1);
   r = cellfun(@length,p1)./diff(eps2,[],2);
    
    rate2(:,:,i) = accumarray([b2 prepost2],r,[],@nanmean,nan);
    
    
     [status,interval] = InIntervals(spks,eps3);
    
   
    
    n1 = histoc(interval(interval>0),1:size(eps3,1));
    
    p1 = mat2cell(spks(status),n1);
   r = cellfun(@length,p1)./diff(eps3,[],2);
    
    rate3(:,:,i) = accumarray([b3 prepost3],r,[],@nanmean,nan);
    
    

end



   %%
   bins = 0 :25:450;
for i = 1:length(uniqueCells)
    
    h = figure('visible','off');
    ax = tight_subplot(1,3);
    
    axes(ax(1))
    spks1 =  spikestimes(spikesIDs(:,3)==i);
    
    [status,interval] = InIntervals(spks1,ineps);
    
    n1 = histoc(interval(interval>0),1:size(ineps,1));

                            
                           %     len_ep = mat2cell(len(status),n1);
    
    
    spks = spks1(status);
    [~,b] = histc(spks,t);
    unit(i).linpos = nan(size(spks,1),1);
    unit(i).linpos(b>0) = len(b(b>0));
    
    
    
    unit(i).pos = mat2cell( unit(i).linpos,n1);
   % set(gca,'ytick',1:25,'yticklabel',1:25)
    % pos = cellfun(@(a) inpaint_nans(a),pos,'uni',0);
    
    
    hold on
     imagesc(bins,[],cell2mat(cellfun(@(a) histoc(a,bins)',pos_inf.len_ep(pos_inf.in_eps),'uni',0))>0)
     colormap([1 1 1 ; .9 .9 .9])
     
   for jj =1:4
         for ii = 1:size(interval_in{jj},1)
             plot([len_in_st{jj}(ii) len_in_fin{jj}(ii)],[interval_in{jj}(ii) interval_in{jj}(ii)],'color',[.7 .7 1],'linewidth',5)
             
         end
     end
    PipeRaster(unit(i).pos)
       
     
         axis on
       title('inbound: stim')
      xlim([0 460])
    axes(ax(2))
     
    [status,interval] = InIntervals(spks1,outeps);
    
n1 = histoc(interval(interval>0),1:size(outeps,1));



    spks = spks1(status);
    [~,b] = histc(spks,t);
unit(i).linpos = nan(size(spks,1),1);
     unit(i).linpos(b>0) = len(b(b>0));

   
          
     unit(i).pos = mat2cell( unit(i).linpos,n1);
     
 %  pos = cellfun(@(a) inpaint_nans(a),pos,'uni',0);

     hold on
     imagesc(bins,[],cell2mat(cellfun(@(a) histoc(a,bins)',pos_inf.len_ep(pos_inf.out_eps),'uni',0))>0)
     colormap([1 1 1 ; .9 .9 .9])
     for jj =1:4
         for ii = 1:size(interval_out{jj},1)
             plot([len_out_st{jj}(ii) len_out_fin{jj}(ii)],[interval_out{jj}(ii) interval_out{jj}(ii)],'color',[.7 .7 1],'linewidth',5)
             
         end
     end
    PipeRaster( unit(i).pos)
     
 xlim([0 460])
     title('outbound: non-stim')
     
      mtit([num2str(uniqueCells(i,:))])
     print(h, '-dpsc2', [basename '.ps'], '-append');
  
  %  uiwait
      close all
end
%%
bins = 5 :15:145;
popmat = nan(length(unit(1).pos),length(bins),length(uniqueCells));
for i = 1:length(uniqueCells)
    
  popmat(:,:,i) =  cell2mat(cellfun(@(a) histoc(a,bins)',unit(i).pos,'uni',0)); 
end
popmat = popmat(:,1:end-1,:);
rate = [];stdev = [];
for i =1:length(uniqueCells)
    rate(:,:,i) = nanmean(linearize(popmat(:,:,i)));
    stdev(:,:,i) = nanstd(linearize(popmat(:,:,i)));
end
popmatz = (popmat - repmat(rate,size(popmat,1),size(popmat,2),1))./repmat(stdev,size(popmat,1),size(popmat,2),1);
for i = 1:size(popmat,2)
    temp(:,:,i) = corrcoef(squeeze(popmatz(:,i,:)),'rows','pairwise');
    
end

figure
imagesc(nanmean(temp,3))
%%

%%

figure

ax = tight_subplot(1,2)
axes(ax(1))
plot(nanmean(rate2(:,1,uniqueCells(:,1) ==2),3),'k')
hold on
plot(nanmean(rate2(:,2,uniqueCells(:,1) ==2),3),'--','color','k')

plot(nanmean(rate2(:,1,uniqueCells(:,1) ==3),3),'r')
plot(nanmean(rate2(:,2,uniqueCells(:,1) ==3),3),'--','color','r')
legend('Sh2 pre','Sh2 post','Sh3 pre','Sh3 post')
title('Shank 2 stim')
xlabel('Stim level')
ylabel('Rate')

axes(ax(2))
plot(nanmean(rate3(:,1,uniqueCells(:,1) ==2),3),'k')
hold on
plot(nanmean(rate3(:,2,uniqueCells(:,1) ==2),3),'--','color','k')

plot(nanmean(rate3(:,1,uniqueCells(:,1) ==3),3),'r')
plot(nanmean(rate3(:,2,uniqueCells(:,1) ==3),3),'--','color','r')
legend('Sh2 pre','Sh2 post','Sh3 pre','Sh3 post')
title('Shank 3 stim')
xlabel('Stim level')
ylabel('Rate')

%%
ineps = cell2mat(cellfun(@(a) [a(1) a(end) ] ,pos_inf.ts_ep(pos_inf.in_eps),'uni',0));
outeps = cell2mat(cellfun(@(a) [a(1) a(end) ] ,pos_inf.ts_ep(pos_inf.out_eps),'uni',0));


t= pos_inf.ts;
kp = InIntervals(sti{2}(:,1),[min(t) max(t)]);
   [status,interval] = InIntervals(t,sti{2}(kp,:));
   len = pos_inf.lin_pos;
   n1 = histoc(interval(interval>0),1:size(sti{2}(kp,:),1));
    len_st = mat2cell( len(status),n1);
    
    
    len_in_st2 = cellfun(@min,len_st(cellfun(@any,len_st)));
    len_in_fin2 = cellfun(@max,len_st(cellfun(@any,len_st)));
       [status,interval_in2] = InIntervals(sti{2}(kp,1),sortrows([ineps;outeps]));
   interval_in2 = interval_in2(interval_in2>0);
   
   kp = InIntervals(sti{4}(:,1),[min(t) max(t)]);
    [status,interval] = InIntervals(t,sti{4}(kp,:));
   len = pos_inf.lin_pos;
   n1 = histoc(interval(interval>0),1:size(sti{4}(kp,:),1));
    len_st = mat2cell( len(status),n1);
    
    
    len_in_st3 = cellfun(@min,len_st(cellfun(@any,len_st)));
    len_in_fin3 = cellfun(@max,len_st(cellfun(@any,len_st)));
       [~,interval_in3] = InIntervals(sti{4}(kp,1),sortrows([ineps;outeps]));
   interval_in3 = interval_in3(interval_in3>0);
   
   %%


close all

figure('positio',[680   142   560   836])
hold on
for i = 1:size(pos_inf.in_eps,1)
    
    if pos_inf.in_eps(i)
       
plot(pos_inf.len_ep{i},pos_inf.ts_ep{i}-pos_inf.ts_ep{i}(1) + i*30,'r','linewidth',2)
    else
        plot(pos_inf.len_ep{i},pos_inf.ts_ep{i}-pos_inf.ts_ep{i}(1) + i*30,'k','linewidth',2)
    end
        
         plot([len_in_st2(interval_in2==i) len_in_fin2(interval_in2==i)],[i*30*ones(sum(interval_in2==i),1) i*30*ones(sum(interval_in2==i),1)],'color',[0 1 1],'linewidth',5)
plot([len_in_st3(interval_in3==i) len_in_fin2(interval_in3==i)],[i*30*ones(sum(interval_in3==i),1) i*30*ones(sum(interval_in3==i),1)],'color',[1 0 1],'linewidth',2)

end

   
set(gca,'yticklabel',[])

xlabel('Position')
ylabel('Trials')
title('Day 1')
%%

st_out = cellfun(@(a) find(a>330 & a <335,1,'first'),pos_inf.len_ep(pos_inf.out_eps),'uni',0);
st_out = cell2mat(cellfun(@(a,b) a(b),pos_inf.ts_ep(pos_inf.out_eps),st_out,'uni',0));
st_in = cellfun(@(a) find(a>330 & a <335,1,'last'),pos_inf.len_ep(pos_inf.in_eps),'uni',0);
st_in = cell2mat(cellfun(@(a,b) a(b),pos_inf.ts_ep(pos_inf.in_eps),st_in,'uni',0));
nChannels =32;
fname_data = 'camkii7_171021.eeg';

st = st_out-5;
dur = 10;
d = nan(length(st),5000);

for i = 1:length(st)
temp =    LoadBinary(fname_data,'nchannels',nChannels,'channels',[ 23],'start',st(i),'duration',dur,'frequency',1250);
    
temp = BandpassFilter(double(temp),1250,[2 300]);
    d(i,:) = temp(3751:end-3750);
end


freqs = logspace(log10(4),log10(300),50);
wt = zeros(50,size(d,2),size(d,1));
%wtCA1_ph = zeros(1250,200,size(d,1),'single');
for j = 1:size(d,1)
   
    [temp,freqlist,psi_array] = awt_freqlist(d(j,:),1250,freqs);
    wt(:,:,j) = (single(abs(temp)'));
 
    
end



[~,interval] = InIntervals(sti{2}(:,1),outeps);
dur = interval(interval>0);


figure
ax = tight_subplot(1,3);
axes(ax(1))
imagesc(-2:1/1250:2,[],nanmedian(wt(:,:,1:dur(1)-1),3))
set(gca,'ydir','normal')
set(gca,'ytick',1:10:50,'yticklabel',round(freqs(1:10:end)))
colormap('jet')
axes(ax(2))
imagesc(-2:1/1250:2,[],nanmedian(wt(:,:,dur),3))
set(gca,'ydir','normal')
set(gca,'ytick',1:10:50,'yticklabel',round(freqs(1:10:end)))

colormap('jet')
axes(ax(3))
imagesc(-2:1/1250:2,[],nanmedian(wt(:,:,dur(end)+1:end),3))
set(gca,'ydir','normal')
set(gca,'ytick',1:10:50,'yticklabel',round(freqs(1:10:end)))
colormap('jet')
%%


st = st_in-5;
dur = 10;
d = nan(length(st),5000);

for i = 1:length(st)
temp =    LoadBinary(fname_data,'nchannels',nChannels,'channels',[ 17],'start',st(i),'duration',dur,'frequency',1250);
    
temp = BandpassFilter(double(temp),1250,[3 300]);
    d(i,:) = temp(3751:end-3750);
end


freqs = logspace(log10(4),log10(300),50);
wt = zeros(50,size(d,2),size(d,1));
%wtCA1_ph = zeros(1250,200,size(d,1),'single');
for j = 1:size(d,1)
   
    [temp,freqlist,psi_array] = awt_freqlist(d(j,:),1250,freqs);
    wt(:,:,j) = single(abs(temp)');
 
    
end



[~,interval] = InIntervals(sti{2}(:,1),outeps);
dur = interval(interval>0);


figure
ax = tight_subplot(1,3);
axes(ax(1))
imagesc(-2:1/1250:2,[],(nanmedian(wt(:,:,1:dur(1)-1),3)))
set(gca,'ydir','normal')
set(gca,'ytick',1:10:50,'yticklabel',round(freqs(1:10:end)))
colormap('jet')
axes(ax(2))
imagesc(-2:1/1250:2,[],nanmedian(wt(:,:,dur),3))
set(gca,'ydir','normal')
set(gca,'ytick',1:10:50,'yticklabel',round(freqs(1:10:end)))
colormap('jet')
axes(ax(3))
imagesc(-2:1/1250:2,[],nanmedian(wt(:,:,dur(end)+1:end),3))
set(gca,'ydir','normal')
set(gca,'ytick',1:10:50,'yticklabel',round(freqs(1:10:end)))
colormap('jet')

