%get position for camkii4 9/12/16

    pix_bin = 140/366;
fname = '/media/sam/SM4/SeqGen/camkii4/20160912/20160912_digitalin.dat';
Fs  =20000;



tt = LoadBinary(fname,'channels',1,'nchannels',1,'frequency',Fs);


ups = find(diff(tt>0)>0);
dwns = find(diff(tt>0)<0);

clear tt




eps = [ups(1)/Fs  ups(end)/Fs];
fs_video = 1./(mode(diff(ups))/Fs);



fname = '/media/sam/SM4/SeqGen/camkii4/20160912/20160912-02.pos';
positions = LoadPositions(fname,30);
positions(positions<=0) = nan;
[t,x,y] = meanCoord(positions);


t(1) = 0;
t = t+eps(1)-19.44;

[x,y] = FixPosition(x,y,t);



%%
[len,pos,xp,yp,del] =linearizePositionTrack(x,y);
 
    len = len*pix_bin;
    
    fils  = getAllExtFiles('/media/sam/SM4/SeqGen/camkii4/20160912/','ait',1);
stims = LoadEvents(fils{1});
sh = nan(length(stims.time/2),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
ok = regexp(stims.description,'\w*','match');
sh  =cellfun(@(a) str2num(a{2}),ok);
amp  =cellfun(@(a) str2num(a{3}),ok);


amp_stim = kmeans(amp(sh==4),7);


%%
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,:);



%%




    %smooth position to get broad epochs of movement in one direction
    sigma = 20;
    kernel  = gaussian2Dfilter([ceil(8*sigma)+1 1],[sigma 1]);
    len1 = nanconvn(len,kernel,'nonanout');
    len = len1;
    
    %find in bound trajectories
    in = diff([0;len1])<0;
    in(end) = 0;
    A = find(diff([0;in])>0);
    B = find(diff([0;in])<0);
    ineps = [A B];
    ep = t(ineps(:,2)) - t(ineps(:,1));
    ineps =  t(ineps(ep>2.5,:)); % must be longer than 3s
    [status,interval] = InIntervals(t,ineps);
    n1 = histoc(interval(interval>0),1:size(ineps,1));
    p1 = mat2cell(len(status),n1);
    kp = cellfun(@range,p1)>50;
    ineps= ineps(kp,:);
    len1(end) = -1;
    %find outbound trajectories
    out = diff([-1;len1])>0;
    
    
        A = find(diff([0;out])>0);
        B = find(diff([0;out])<0);
    outeps = [A B];
    
       ep = t(outeps(:,2)) - t(outeps(:,1));
    outeps =  t(outeps(ep>2.5,:)); % must be longer than 3s
     [status,interval] = InIntervals(t,outeps);
    n1 = histoc(interval(interval>0),1:size(outeps,1));
    p1 = mat2cell(len(status),n1);
    kp = cellfun(@range,p1)>50;
    outeps= outeps(kp,:);
    

    
    
    
   [status,interval] = InIntervals(t,ineps);
   n1 = histoc(interval(interval>0),1:size(ineps,1));
    len1 = mat2cell( len(status),n1);
    
       [status,interval] = InIntervals(t,outeps);
   n1 = histoc(interval(interval>0),1:size(outeps,1));
    len2 = mat2cell( len(status),n1);
    
    
    
    
    
     [status,interval] = InIntervals(t,sti{3});
   n1 = histoc(interval(interval>0),1:size(sti{3},1));
    len_st = mat2cell( len(status),n1);
    
    
    len_in_st = cellfun(@min,len_st);
    len_in_fin = cellfun(@max,len_st);
       [~,interval_in] = InIntervals(sti{3}(:,1),ineps);
   interval_in = interval_in(interval_in>0);
    
%%

bins = 0 :5:160;

             
for i = 1:length(uniqueCells)
    
    h = figure('visible','off');
    ax = tight_subplot(1,2);
    
    axes(ax(1))
   spks1 =  spikestimes(spikesIDs(:,3)==i);
    
    [status,interval] = InIntervals(spks1,ineps);
    
n1 = histoc(interval(interval>0),1:size(ineps,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)',len1,'uni',0))>0)
     colormap([1 1 1 ; .9 .9 .9])
     
     for ii = 1:size(interval_in,1)
          plot([len_in_st(ii) len_in_fin(ii)],[interval_in(ii) interval_in(ii)],'color',[.7 .7 1],'linewidth',5)
          
     end
    PipeRaster(unit(i).pos)
       
 %PipeRaster(len1)
   
     
      
       title('inbound: stim')
      xlim([0 160])
    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));
     [~,b] = bestmatch(ineps(:,1),outeps(:,1));
       [~,interval] = InIntervals(sti{4}(:,1),ineps);
     interval_out = cell2mat(cellfun(@(a) find(outeps(:,1)>a,1) ,num2cell(ineps(interval(interval>0)),2),'uni',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)',len2,'uni',0))>0)
     colormap([1 1 1 ; .9 .9 .9])
       
          
            for ii = 1:size(interval_out,1)
          plot([len_in_st(ii) len_in_fin(ii)],[interval_out(ii) interval_out(ii)],'color',[.7 .7 1],'linewidth',5)
          
            end
     
    PipeRaster( unit(i).pos)
     
 xlim([0 160])
     title('outbound: non-stim')
      
      mtit([num2str(uniqueCells(i,:))])
   %  print(h, '-dpsc2', 'camkii4_91116.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))
%%



pre = sti{4}(sti{4}(:,2)<sti{3}(1,1) ,:);

co = sti{3};

post = sti{4}(sti{4}(:,2)>sti{3}(end,1),:);

bins = -1:.01:1;
rate_pre = nan(size(uniqueCells,1),length(bins));

for i = 1:size(uniqueCells,1)
    spks = spikestimes(spikesIDs(:,3)==i);
    [status,interval] = InIntervals(spks,[pre(:,1)-1 pre(:,2)+1]);
    
   
    
    n1 = histoc(interval(interval>0),1:size(pre,1));
    
    p1 = mat2cell(spks(status),n1);

    p1 = cellfun(@(a,b) a-b,p1,num2cell(pre(:,1),2),'uni',0);
   rate_pre(i,:) =  histc(cell2mat(p1),bins)/length(pre)/.01;
end

rate_co = nan(size(uniqueCells,1),length(bins));

for i = 1:size(uniqueCells,1)
    spks = spikestimes(spikesIDs(:,3)==i);
    [status,interval] = InIntervals(spks,[co(:,1)-1 co(:,2)+1]);
    
   
    
    n1 = histoc(interval(interval>0),1:size(co,1));
    
    p1 = mat2cell(spks(status),n1);

    p1 = cellfun(@(a,b) a-b,p1,num2cell(co(:,1),2),'uni',0);
   rate_co(i,:) =  histc(cell2mat(p1),bins)/length(co)/.01;
end

rate_post = nan(size(uniqueCells,1),length(bins));

for i = 1:size(uniqueCells,1)
    spks = spikestimes(spikesIDs(:,3)==i);
    [status,interval] = InIntervals(spks,[post(:,1)-1 post(:,2)+1]);
    
   
    
    n1 = histoc(interval(interval>0),1:size(post,1));
    
    p1 = mat2cell(spks(status),n1);

    p1 = cellfun(@(a,b) a-b,p1,num2cell(post(:,1),2),'uni',0);
   rate_post(i,:) =  histc(cell2mat(p1),bins)/length(post)/.01;
end

%%

sh = repmat([find(diff([uniqueCells(:,1)]))],1,2)+1;
ax = tight_subplot(1,3);

ok = [rate_pre rate_co rate_post];
ok = zscore(ok,[],2);

rate_prez = ok(:,1:size(rate_pre,2));

rate_coz = ok(:,size(rate_pre,2)+1:2*size(rate_pre,2));

rate_postz =  ok(:,2*size(rate_pre,2)+1:end);


axes(ax(1))



imagesc(bins,[],nanzscore(rate_pre,[],2),[0 10])
hold on
plot(repmat([bins(1) bins(end)],2,1)',sh'-.5,'w')
title('Pre (sh3)')
axes(ax(2))


imagesc(bins,[],nanzscore(rate_co,[],2),[0 10])
hold on
plot(repmat([bins(1) bins(end)],2,1)',sh'-.5,'w')

title('Tack Pairing (n=13, sh3/4)')
axes(ax(3))

imagesc(bins,[],nanzscore(rate_post,[],2),[0 10])



hold on
plot(repmat([bins(1) bins(end)],2,1)',sh'-.5,'w')

title('Post (sh3)')
