%% CANNOT LOOK AT SPECTRUM SINCE THE STIM LEVEL CHANGES FOR EVERY PULSE

% JC.name = '/media/sam/SM4/SeqGen/JS13/114618_150118_150343';
% JC.JCID = [1 15];
% JC.post(1).ID = [1 13];
% JC.post(2).ID = [1 14];
%

cd /media/sam/SM4/SeqGen/JS13/20161005_161005_133644/
clear all
% get ripples
dirName = pwd;
fils  = getAllExtFiles(dirName,'rip',1);

ripples = LoadEvents(fils{1});

rip_on =  ripples.time(cellfun(@(a) any(regexp(a,'Ripple start')), ...
    ripples.description));

rip_off =  ripples.time(cellfun(@(a) any(regexp(a,'Ripple stop')), ...
    ripples.description));

rip_eps = [rip_on rip_off];



SetCurrentSession('/media/sam/SM4/SeqGen/JS13/20161005_161005_133644/20161005_161005_133644.xml')
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,:);


jxta = [1 15];

intID = [1 13;1 14];
badPYR = [1 3;1 6];


MUA = setdiff(units,[1 1;jxta;intID;badPYR],'rows');
pre_spks = spikestimes(ismember(spikesIDs(:,1:2),jxta,'rows'));


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

for j = 1
    
    sti= 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(dirName,'aip',1);
[amp,sh] = getPulseAmplitude(fils{1});

PETH = getPETH(pre_spks,sti(:,1),.2,.2,.001);

kp = diff(sti,[],2)>.01;

JC.PETH = PETH(kp,:);
JC.all_sti = sti(kp,:);




kp = mean(PETH(:,200:300),2)>5;
JC.good_sti = sti(kp,:);
epochs = [sti(kp,1) sti(kp,2)];
%epochs = excludeEpochs(epochs,rip_eps);
adjustedSpikes = cellfun(@(a) pre_spks(pre_spks>=a(1) & pre_spks<= a(2)),num2cell(epochs,2),'uni',0);
st = epochs(cellfun(@any,adjustedSpikes),1);
adjustedSpikes = adjustedSpikes(cellfun(@any,adjustedSpikes));
earlySpikes = cell2mat(cellfun(@(a,b) a((a-b)<.05),adjustedSpikes,num2cell(st,2),'uni',0));
lateSpikes = cell2mat(cellfun(@(a,b) a((a-b)>.05),adjustedSpikes,num2cell(st,2),'uni',0));
firstSpikes = cellfun(@(a) a(1),adjustedSpikes);
timFirstSpikes = cellfun(@(a,b) a(1)-b,adjustedSpikes,num2cell(st,2));
adjustedSpikes = cell2mat(adjustedSpikes);
JC.firstSpikeTiming = firstSpikes-st;






[status,interval] = InIntervals(pre_spks,sti);

[status1,interval] = InIntervals(pre_spks,rip_eps);


spont_spike = pre_spks(~status& ~status1);



for ii = 1:size(intID,1)+1
    
    if ii<=size(intID,1)
        post_spks = spikestimes(ismember(spikesIDs(:,1:2),intID(ii,:),'rows'));
    else
        post_spks = spikestimes(ismember(spikesIDs(:,1:2),MUA,'rows'));
    end
    
    
    for kk = 1:5
        
        
        
        switch kk
            case 1
                pre = spont_spike;
                field = 'spont';
                description = 'CCG for spontanous spikes';
            case 2
                pre = adjustedSpikes;
                field = 'evoked';
                description = 'all evoked spikes for effective stimulations';
            case 3
                pre = firstSpikes;
                field = 'first';
                description = 'all first evoked spikes for effective stimulations';
            case 4
                pre = earlySpikes;
                field = 'early';
                description = 'all evoked spikes for effective stimulations before 50ms';
            case 5
                pre = lateSpikes;
                field = 'late';
                description = 'all evoked spikes for effective stimulations after 50ms';
                
        end
        
        
        spk_tims_1 = [pre;post_spks];
        spk_IDs_1 = [ones(length(pre),1);2*ones(length(post_spks),1)];
        [spk_tims_1,b] = sort(spk_tims_1);
        spk_IDs_1 = spk_IDs_1(b);
        temp1 = CCG(spk_tims_1, spk_IDs_1,'binsize',.0008,'duration',.2);
        
        [~,pred] = cch_conv(temp1(:,1,2),round(.015/.0008));
        
        
        prob_uncor = temp1(:,1,2)/length(pre);
        
        if ii<=size(intID,1)
            JC.post(ii).(field).raw = temp1(:,1,2);
            JC.post(ii).(field).prob_uncor = prob_uncor;
            JC.post(ii).(field).prob = (temp1(:,1,2)-pred)/length(pre);
            JC.post(ii).(field).n = length(pre);
            JC.post(ii).(field).description = description;
            JC.post(ii).(field).trans = nansum(JC.post(ii).(field).prob(127:129));
        else
            
            JC.MUA.(field).raw = temp1(:,1,2);
            JC.MUA.(field).prob_uncor = prob_uncor;
            JC.MUA.(field).prob = (temp1(:,1,2)-pred)/length(pre);
            JC.MUA.(field).n = length(pre);
            
            JC.MUA.(field).trans = nansum(JC.MUA.(field).prob(127:129));
            
            
        end
        
    end
end




p1 = [JC.post(1).spont.trans JC.post(1).evoked.trans JC.post(1).first.trans JC.post(1).early.trans JC.post(1).late.trans];
p2 = [JC.post(2).spont.trans JC.post(2).evoked.trans JC.post(2).first.trans JC.post(2).early.trans JC.post(2).late.trans];
save('JC_CCG.mat','JC')

%%


%%
figure
%first spike
epochs = [sti(kp,1) sti(kp,1)+.025];
adjustedSpikes = cellfun(@(a) pre_spks(pre_spks>=a(1) & pre_spks<= a(2))-(a(1)),num2cell(epochs,2),'uni',0);
adjustedSpikes = adjustedSpikes(cellfun(@any,adjustedSpikes));

PipeRaster(adjustedSpikes)
%%

kp = diff(sti,[],2)>.05 & amp>3500 & sti(:,1)<12911 & sti(:,1)>12741;
for i = max(spikesIDs(:,3))
    figure ;temp = spikestimes(spikesIDs(:,3) == i);
    PETH = getPETH(temp,sti(kp,1),.12,.19,.001);
    bar(-.12:.001:.19,nanmean(PETH))
    title(num2str(uniqueCells(i,:)))
    
end

kp = diff(sti,[],2)>.05 & amp>3500 & sti(:,1)<12911 & sti(:,1)>12741;
epochs = [sti(kp,1)-.12 sti(kp,1)+.19];

adjustedSpikes = cellfun(@(a) pre_spks(pre_spks>=a(1) & pre_spks<= a(2))-(a(1)+.12),num2cell(epochs,2),'uni',0);



epochs = [sti(2000:3000,1)-.1 sti(2000:3000,1)+.1];

adjustedSpikes = cellfun(@(a) pre_spks(pre_spks>=a(1) & pre_spks<= a(2))-(a(1)+.1),num2cell(epochs,2),'uni',0);
figure
PipeRaster(adjustedSpikes)


