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

clear all
cd /media/sam/SM4/SeqGen/RM14/20170711g/
SetCurrentSession('/media/sam/SM4/SeqGen/RM14/20170711g/20170711g.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 33];

intID = [1 27;1 28; 1 29; 1 30; 1 31; 1 32];
badPYR = [ ];


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

dirName = pwd;
rip_eps = [0 .01];
% get ripples

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


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

for j = 7
    
    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.all_sti = sti(kp,:);
JC.PETH = PETH(kp,:);
kp = mean(PETH(:,200:250),2)>5 & diff(sti,[],2)>.01;
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(128:130));
        
         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


save('JC_CCG.mat','JC')

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


 ses.ID = spikesIDs(:,1:2);
ses.juxta_ID = jxta;
 ses.times = spikestimes;
 ses.stimTimes = sti;
 ses.name = pwd;
 ses.rippleTimes = rip_eps;
 save('JC_groundTruth.mat','ses')


%p3 = [JC.post(3).spont.trans JC.post(3).evoked.trans JC.post(3).first.trans JC.post(3).early.trans JC.post(3).late.trans];
%p4= [JC.post(4).spont.trans JC.post(4).evoked.trans JC.post(4).first.trans JC.post(4).early.trans JC.post(4).late.trans];

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

kp = mean(PETH(:,200:250),2)>5 & diff(sti,[],2)>.01;
epochs = [sti(kp,1)-.1 sti(kp,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)


%%


kp = diff(sti,[],2)>.015 ;
epochs = [sti(kp,1)-.1 sti(kp,1)+.1];

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


%%

%get post triggered lfp

d = LoadBinary('20170615c_EC.dat','nchannels',33,'channels',1,'frequency',30000);
%%
[status,interval] = InIntervals(post_spks,[sti(:,1)-.5 sti(:,2)+.5]);
spk_ts = post_spks(~status);
spk_ts = spk_ts(diff([0;spk_ts])>.05);
Fs = 30000;
before = .2;
after = .2;
before = round(before*Fs);
after = round(after*Fs);
dur = after+before+1;
spk_ts = round(spk_ts*Fs);
idx = repmat(-before:after,size(spk_ts,1),1) + repmat(spk_ts,1,dur);
kp = all(idx>0,2) & all(idx<=length(d),2);
idx = idx(kp,:);
lfp = d(idx);
%%

k = designfilt('bandstopiir','FilterOrder',2, ...
    'HalfPowerFrequency1',59,'HalfPowerFrequency2',61, ...
    'DesignMethod','butter','SampleRate',Fs);
for i = 1:size(lfp)
    lfp1(i,:) = filtfilt(k,double(lfp(i,:)));
end

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

adjustedSpikes = cellfun(@(a) a(1),adjustedSpikes);
bins = 0:.0005:.025;
bar(bins, histc(adjustedSpikes,bins))
%%


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

spont_spike = pre_spks(~status);

spk_tims_1 = [spont_spike;post_spks];
spk_IDs_1 = [ones(length(spont_spike),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);
spont_CCG = temp1(:,1,2)/length(spont_spike);

%%


epochs = [sti(kp,1) sti(kp,2)];
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));
ealySpikes = cell2mat(cellfun(@(a,b) a((a-b)<.015),adjustedSpikes,num2cell(st,2),'uni',0));
lateSpikes = cell2mat(cellfun(@(a,b) a((a-b)>.015),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);
st = firstSpikes-st;
%%


spk_tims_1 = [adjustedSpikes;post_spks];
spk_IDs_1 = [ones(length(adjustedSpikes),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);
first_CCG = temp1(:,1,2)/length(adjustedSpikes);
plot(-.1:.0008:.1,first_CCG)
hold on
plot([0 0],[0 .02])
xlim([-.05 .05])
plot(-.1:.0008:.1,spont_CCG,'k'    )

%%


latefirst = firstSpikes(timFirstSpikes>.01);

temp = [latefirst;post_spks];
tempid = [ones(length(latefirst),1);2*ones(length(post_spks),1)];

[temp,b] = sort(temp);
tempid = tempid(b);


dt = .02;idx = 1;ev =[];
%loop over spikes
for i = 1:length(temp)
    
    
    %count back
    j = -1;
    
    while i+j > 0 &&  tempid(i) ==1 && (temp(i) - temp(i+j)) <= dt && tempid(i+j) ==2
        
        ev{idx} = temp(i+j);
        idx = idx+1;
        j = j-1;
        
    end
    
    
end



%%


% plot spike triggered spectra

%load data
d = LoadBinary('114618_150118_150343_EC.dat','nchannels',16,'channels',9,'frequency',20000);

%%
ok = find(kp);
%get artefact
on = round(sti(ok(1:3000),1)*20000);

idx = repmat(-2000:10000,size(on,1),1) + repmat(on(:,1),1,12001);
dat = d(idx);

d = LoadBinary('114618_150118_150343_EC.dat','nchannels',16,'channels',9,'frequency',20000);


a =  LoadBinary('114618_150118_150343_analogin.dat','nchannels',2,'channels',0,'frequency',20000);
art = a(idx);
%subtract the mean artefact from the lfp
artefact = int16(nanmedian(dat)');
%%
for i = 1:length(on)
    
    d(idx(i,:)) = d(idx(i,:)) - artefact;
end
%%

f = logspace(log10(3),log10(300),50);
spks1 = spikes(:,1);
[status,interval] = InIntervals(spks1,sti);
pre_spks = spks1(status);
% spks = spks1(~(spks1>sti(1,1) & spks1<sti(end,2)));
% spks = spks(randsample(1:length(spks),size(lfp1,1)));
S_stim = GetSpikeTriggeredSpectra(pre_spks,d,3,3,1250,f);
%%
figure
imagesc(-2:1/1250:2,[],nanzscore(nanmean(S_stim(1250:end-1250,:,:),3))',[-1 1])
set(gca,'ydir','normal')
colormap('jet')
set(gca,'ytick',1:10:50,'yticklabel',round(f(1:10:end)))
hold on
plot([0 0],[0 50],'w')
plot([.015 .015],[0 50],'w')

%%
for i = 1:1000
    S_stim_cont(i).u = nanmean(GetSpikeTriggeredSpectra(pre_spks+rand(size(pre_spks,1),1)*60,d,3,3,1250,f),3);
    S_stim_cont(i).v = nanvar(GetSpikeTriggeredSpectra(pre_spks+rand(size(pre_spks,1),1)*60,d,3,3,1250,f),[],3);
    i
end
B1 = cell2mat(cellfun(@(x) permute(x,[3 1 2]),{S_stim_cont.u},'UniformOutput',false)');
%%
p = nan(size(S_stim_cont,1)-2500,size(S_stim_cont,2));
for i = 1:size(S_stim_cont,1)-2500
    for j = 1:size(S_stim_cont,2)
        [h,p(i,j)] =  ttest2(log10(squeeze(S_stim_cont(1250+i,j,:))),log10(squeeze(S_stim(1250+i,j,:))));
    end
end


%%
figure
imagesc( -2:1/1250:2,[],-log10(p)')
colormap('jet')
set(gca,'ydir','normal')
colormap('jet')
set(gca,'ytick',1:10:50,'yticklabel',round(f(1:10:end)))
hold on
plot([0 0],[0 50],'w')
plot([.015 .015],[0 50],'w')


%%
spks1 = spikes(:,1);
[status,interval] = InIntervals(spks1,[sti(:,1)-.5 sti(:,end)+.5]);
pre_spks = spks1(status);

pre_spks = pre_spks(randsample(1:length(pre_spks),size(S_stim,3)));
S_nostim = GetSpikeTriggeredSpectra(pre_spks,d,3,3,1250,f);


%%


S_nostim_cont = GetSpikeTriggeredSpectra(pre_spks+rand(size(pre_spks))*50,d,3,3,1250,f);

%%
p_nostim = nan(size(S_nostim_cont,1)-2500,size(S_nostim_cont,2));
for i = 1:size(S_nostim_cont,1)-2500
    for j = 1:size(S_nostim_cont,2)
        p_nostim(i,j) =  ranksum(squeeze(S_nostim_cont(1250+i,j,:)),squeeze(S_nostim(1250+i,j,:)));
    end
end


%%
figure
imagesc( -2:1/1250:2,[],-log10(p_nostim)',[0 5])
colormap('jet')
set(gca,'ydir','normal')
colormap('jet')
set(gca,'ytick',1:10:50,'yticklabel',round(f(1:10:end)))
xlim([-1 1])
%%
figure
imagesc( -2:1/1250:2,[],-log10(p)',[0 5])
colormap('jet')
set(gca,'ydir','normal')
colormap('jet')
set(gca,'ytick',1:10:50,'yticklabel',round(f(1:10:end)))
xlim([-1 1])

%%


ok = nanzscore(cat(3,S_stim(1250:end-1250,:,:),S_nostim(1250:end-1250,:,:)));
S_stim1 = ok(:,:,1:size(S_stim,3));
S_nostim1 = ok(:,:,size(S_stim,3)+1:end);


%%
figure
imagesc(-2:1/1250:2,[],nanmean(S_stim_cont(1250:end-1250,:,:),3)',[-.1 .3])
set(gca,'ydir','normal')
colormap('jet')
set(gca,'ytick',1:10:50,'yticklabel',round(f(1:10:end)))
hold on
plot([0 0],[0 50],'w')
plot([.015 .015],[0 50],'w')

figure
imagesc(-2:1/1250:2,[],nanmean(S_nostim1(1250:end-1250,:,:),3)' ,[-.1 .3])
set(gca,'ydir','normal')
colormap('jet')
set(gca,'ytick',1:10:50,'yticklabel',round(f(1:10:end)))

set(gca,'ytick',1:10:50,'yticklabel',round(f(1:10:end)))


%%



%%


figure
imagesc(-2:1/1250:2,[],nanmean(spect_filt_spk1,3)' ,[-.2 .4])
set(gca,'ydir','normal')
colormap('jet')
xlim([-1 1])
hold on
plot([0 0],[0 50],'w')
plot([.015 .015],[0 50],'w')
set(gca,'ytick',1:10:50,'yticklabel',round(f(1:10:end)))

figure
imagesc(-2:1/1250:2,[],nanmean(spect_nostim_spk1,3)',[-.2 .4])
set(gca,'ydir','normal')
colormap('jet')
set(gca,'ytick',1:10:50,'yticklabel',round(f(1:10:end)))


%%
f = logspace(log10(3),log10(300),50);
spect_raw = nan(6001,50,size(dat,1));
spect_filt = nan(6001,50,size(dat,1));
for i = 1:size(dat,1)
    [temp] = awt_freqlist(double(dat(i,1:10:end)),3000,f);
    
    spect_raw(:,:,i) = abs(temp);
    [temp] = awt_freqlist(double(dat1(i,1:10:end)),3000,f);
    spect_filt(:,:,i) = abs(temp);
end


figure
imagesc(nanzscore(nanmean(spect_raw,3)',[],2),[0 2])
set(gca,'ydir','normal')
colormap('jet')
%%
figure
imagesc(nanzscore(nanmean(spect_filt,3)',[],2),[0 2])
set(gca,'ydir','normal')
colormap('jet')
%%
data1 = double(dat1(:,1:10:end))';
data2 = getPETH(spikes(:,1),sti(:,1),1,1,1/3000)';
params.Fs = 3000;
params.fpass = [4 300];
params.tapers = [3 5];
[C,phi,S12,S1,S2,t,f,zerosp]=cohgramcpb(data2(:),data2(:),[.5 1/300],params);


