

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);


%%

dirName = pwd;
basename = '20161201';
%find CCGs for different stim levels


 mono_res = MonoSynConvClick3 (spikesIDs,spikestimes,'plot',false);
    nCel = size(mono_res.completeIndex,1);
    uniqueCells = mono_res.completeIndex(:,1:2);
  
    ID = ones(26,1);
    ID([21 24 25]) = 4;
   
    pre_post = mono_res.sig_con(all([mono_res.n(mono_res.sig_con(:,1)) mono_res.n(mono_res.sig_con(:,2))] >500,2),:);
  %%
    CCGxTime = CCG_stim_triplet2(dirName,basename, pre_post,ID,{[],[st1(:,1) st1(:,1)+.05]},units);

%make histograms


%%
figure
ts = logspace(log10(.5),log10(2000),40)/1000;
ts = ts(1:end-1);
ses(1).CCGxtime_triplet =CCGxTime;


%unpack triplet
trans = [];
pre1 = [];
pre2 = [];
%ID = [];
stf = [];
n=[];
same_pre_cell = [];
same_post_cell = [];
common_instar = [];
common_outstar = [];
second_instar = [];postID=[];
ID = [];
celID1 = [];celID2 =[];
ccgR =[];
for i = 1:length(ses)
    
    
   
    for j = 1:length(ses(i).CCGxtime_triplet.spont.post)
        
         
                    
        for l = 1:length(ses(i).CCGxtime_triplet.spont.post{j}.pre)
            
            if ~isempty(ses(i).CCGxtime_triplet.spont.post{j}.pre(l).triplet) && ...
                isfield(ses(i).CCGxtime_triplet.spont.post{j}.pre(l).triplet(1),'prob')
                
            ID_cel = cell2mat({ses(i).CCGxtime_triplet.spont.post{1}.pre.ID});
                for k = 1:length(ses(i).CCGxtime_triplet.spont.post{j}.pre(l).triplet)
                    
               
                    
                    
                    numspk = ses(i).CCGxtime_triplet.spont.post{j}.pre(l).triplet(k).numspk;
                    if numspk(1)>0
                        trans = [trans;nansum(ses(i).CCGxtime_triplet.spont.post{j}.pre(l).triplet(k).prob(201:203,:))];
                        pre1=  [pre1;nansum(ses(i).CCGxtime_triplet.spont.post{j}.pre(l).prob(201:203,:))];
                  
                         pre_ID = ses(i).CCGxtime_triplet.spont.post{j}.pre(l).triplet(k).name(3);
                          pID = ses(i).CCGxtime_triplet.spont.post{j}.pre(pre_ID).ID;
                       %
                        same_pre_cell = [same_pre_cell;ses(i).CCGxtime_triplet.spont.post{j}.pre(l).triplet(k).same_pre_cell];
                        same_post_cell = [same_post_cell;ses(i).CCGxtime_triplet.spont.post{j}.pre(l).triplet(k).same_post_cell];
                        common_instar = [common_instar;ses(i).CCGxtime_triplet.spont.post{j}.pre(l).triplet(k).common_instar];
                
                         pre2=  [pre2;ses(i).CCGxtime_triplet.spont.post{j}.pre(pre_ID).trans];
                      
                        common_outstar = [common_outstar;ses(i).CCGxtime_triplet.spont.post{j}.pre(l).triplet(k).common_outstar];
                        second_instar = [second_instar;ses(i).CCGxtime_triplet.spont.post{j}.pre(l).triplet(k).second_instar];
                     
                      postID = [postID;ID_cel(ses(i).CCGxtime_triplet.spont.post{j}.name(3))];
                        ID=  [ID;pID];
                          n = [n;numspk'];
                        ccgR = [ccgR;ses(i).CCGxtime_triplet.spont.post{j}.pre(l).raw'];
                        celID1 = [celID1;[i ses(i).CCGxtime_triplet.spont.post{j}.pre(l).name]];
                         celID2 = [celID2;[i ses(i).CCGxtime_triplet.spont.post{j}.pre(l).triplet(k).name]];
                    end
                end
            end
        end
    end
end

trans(n<10) = nan;

temp1 = (trans(same_post_cell==1,:)-pre1(same_post_cell==1,:))./pre1(same_post_cell==1,:);

%unpack triplet
trans = [];
pre1 = [];
pre2 = [];
%ID = [];
stf = [];

same_pre_cell = [];
same_post_cell = [];
common_instar = [];
common_outstar = [];
second_instar = [];postID=[];
ID = [];
celID1 = [];celID2 =[];
ccgR =[];n=[];
for i = 1:length(ses)
    
    
   
    for j = 1:length(ses(i).CCGxtime_triplet.evoked.post)
        
         
                    
        for l = 1:length(ses(i).CCGxtime_triplet.evoked.post{j}.pre)
            
            if ~isempty(ses(i).CCGxtime_triplet.evoked.post{j}.pre(l).triplet) && ...
                isfield(ses(i).CCGxtime_triplet.evoked.post{j}.pre(l).triplet(1),'prob')
                
            ID_cel = cell2mat({ses(i).CCGxtime_triplet.evoked.post{1}.pre.ID});
                for k = 1:length(ses(i).CCGxtime_triplet.evoked.post{j}.pre(l).triplet)
                    
               
                    
                    
                    numspk = ses(i).CCGxtime_triplet.evoked.post{j}.pre(l).triplet(k).numspk;
                    if numspk(1)>0
                        trans = [trans;nansum(ses(i).CCGxtime_triplet.evoked.post{j}.pre(l).triplet(k).prob(201:203,:))];
                        pre1=  [pre1;nansum(ses(i).CCGxtime_triplet.evoked.post{j}.pre(l).prob(201:203,:))];
                  
                         pre_ID = ses(i).CCGxtime_triplet.evoked.post{j}.pre(l).triplet(k).name(3);
                          pID = ses(i).CCGxtime_triplet.evoked.post{j}.pre(pre_ID).ID;
                       %
                        same_pre_cell = [same_pre_cell;ses(i).CCGxtime_triplet.evoked.post{j}.pre(l).triplet(k).same_pre_cell];
                        same_post_cell = [same_post_cell;ses(i).CCGxtime_triplet.evoked.post{j}.pre(l).triplet(k).same_post_cell];
                        common_instar = [common_instar;ses(i).CCGxtime_triplet.evoked.post{j}.pre(l).triplet(k).common_instar];
                
                         pre2=  [pre2;ses(i).CCGxtime_triplet.evoked.post{j}.pre(pre_ID).trans];
                      n = [n;numspk'];
                        common_outstar = [common_outstar;ses(i).CCGxtime_triplet.evoked.post{j}.pre(l).triplet(k).common_outstar];
                        second_instar = [second_instar;ses(i).CCGxtime_triplet.evoked.post{j}.pre(l).triplet(k).second_instar];
                     
                      postID = [postID;ID_cel(ses(i).CCGxtime_triplet.evoked.post{j}.name(3))];
                        ID=  [ID;pID];
                        ccgR = [ccgR;ses(i).CCGxtime_triplet.evoked.post{j}.pre(l).raw'];
                        celID1 = [celID1;[i ses(i).CCGxtime_triplet.evoked.post{j}.pre(l).name]];
                         celID2 = [celID2;[i ses(i).CCGxtime_triplet.evoked.post{j}.pre(l).triplet(k).name]];
                    end
                end
            end
        end
    end
end

trans(n<10) = nan;


temp = (trans(same_post_cell==1,:)-pre1(same_post_cell==1,:))./pre1(same_post_cell==1,:);
semilogx(ts(ts>.004 & ts<.1),temp(:,ts>.004& ts<.1),'o','color','k')
hold on
semilogx(ts(ts>.004),nanmedian(temp1(:,ts>.004)),'r','linewidth',3)
semilogx([.003 2],[0 0],'k')


%%
figure
hold on
for i = 1:3
    
eps = [st1(cat1==i,1)  st1(cat1==i,1) + .05];

[status,interval] = InIntervals(spikestimes,eps);


spikestimes_kp = spikestimes(status);
spikesIDs_kp = spikesIDs(status,3);
[ccgR1] = CCG(spikestimes_kp,spikesIDs_kp,'binSize',.0004,'duration',.1);
ccg =[];
for j =  1:size(pre_post,1)
    n = sum(spikesIDs_kp ==pre_post(j,1));
   ccg = [ccg;ccgR1(:,pre_post(j,1),pre_post(j,2))'/n];
end
plot(nanmean(ccg))
end