batchVettedRestricted = StartBatch(@BatchGetDeltasPostSleep,'nonstim.batch');
X = get(batchVettedRestricted,'UserData');

%% Compute correlations between ripple activity and following corical activity (200ms window)

movingWindow = -0.5:0.02:0.5;
t = movingWindow;
nSessions = size(X,1);
duration = 0.2;
nShifts = 5;

[mc0, mc1, mp0, mp1, mi0, mi1] = deal([]);

for j=1:2,
	for ses=1:nSessions,
		[deltas,occ,groups,pfcUnit,ripples,pfc,hpc,~,session] = X{ses,:};
		[sessionPath,sessionName] = fileparts(session);
		occ = occ>0.5;
		groups = round(groups); if min(groups)>8,groups=groups-8;end
		peaks = X{ses,7};
		pfcUnit = pfcUnit(:,1); if min(pfcUnit)>8,pfcUnit=pfcUnit-8;end
		
		% 		ripples = dlmread(fullfile(sessionPath,'ripples'));
		ripples = dlmread(['/home/raly/projects/delta/RespondingToReviewers/NewRipples/' sessionName '-allRipples']);
		
		positions = dlmread(fullfile(fileparts(session),'positions'));
		behaviour = positions([1 end],1)';
		ripples = Restrict(ripples,[behaviour(end) Inf]);
		preceding = InIntervals(ripples(:,2),[deltas(:,2)-0.2 deltas(:,2)]);
		
		switch j
			case 1,
				ok = ~preceding;
			case 2
				ok = preceding;
		end
		
		if sum(ok)<20 | sum(~ok)<20,
			nCorrs(ses,j) = 0;nEvents(ses,j) = 0;continue;
		end
		
		ripples = ripples(ok,:);
		rIntervals = [ripples(:,2)-0.1 ripples(:,2)+0.1]; % special: make them 200 ms for a change
		source = CountSpikesIntervals(hpc,rIntervals);
		cc = []; pp = [];
		for i=1:length(movingWindow),
			intervals = bsxfun(@plus,ripples(:,2),[-1 1]*duration/2+movingWindow(i));
			target = CountSpikesIntervals(pfc,intervals);
			thisc = []; thisp = [];
			for shift = 0:nShifts,
				[thiscShifted,thispShifted] = nancorr(target((1:end-nShifts)+shift,:),source((1:end-nShifts),:),'type','spearman'); thisc = [thisc thiscShifted(:)]; thisp = [thisp thispShifted(:)];
			end
			cc = cat(3,cc,thisc); pp = cat(3,pp,thisp);
		end
		[idSource,idTarget] = meshgrid(1:size(source,2),1:size(target,2));
		switch j
			case 1,
				mc0 = cat(1,mc0,cc); mp0 = cat(1,mp0,pp); mi0 = cat(1,mi0,[idSource(:) idTarget(:)]);
			case 2
				mc1 = cat(1,mc1,cc); mp1 = cat(1,mp1,pp); mi1 = cat(1,mi1,[idSource(:) idTarget(:)]);
		end
		nCorrs(ses,j) = numel(idSource);
		nEvents(ses,j) = size(source,1);
		ses
	end
end

%% Plot comparison 

figure(1)

idx = FindClosest(t,0.13);

subplot(2,3,1);
g = Group(mc0(:,1,idx),mc0(:,2,idx));
anovabar(g(:,1),g(:,2));
set(gca,'xtick',1:2,'xticklabel',{'data (no delta)','shifted control'});
ylabel('mean corr between [-100 100]ms ripple window and [0 200]ms cortical spikes');
title(['No delta, p=' num2str(signrank(g(g(:,2)==1),g(g(:,2)==2)))]);

subplot(2,3,2);
g = Group(mc1(:,1,idx),mc1(:,2,idx));
anovabar(g(:,1),g(:,2));
set(gca,'xtick',1:2,'xticklabel',{'data (delta)','shifted control'});
title(['Delta, p=' num2str(ranksum(g(g(:,2)==1),g(g(:,2)==2)))]);

subplot(2,3,3);
g = Group(mc0(:,1,idx),mc1(:,1,idx));
anovabar(g(:,1),g(:,2));
title(['No delta vs delta, p=' num2str(signrank(g(g(:,2)==1),g(g(:,2)==2)))]);
set(gca,'xtick',1:2,'xticklabel',{'no delta','delta'});

EquateScales(1,2,3);

% clear sig
for j=1:2,
	switch j
		case 1, mc = mc0; color = 'b'; titl = 'NoDelta';
		case 2, mc = mc1; color = 'r';  titl = 'delta within 200 ms';
	end
	sesh = repelem((1:nSessions)',nCorrs(:,j));
	subplot(2,2,2+j+(i-1)*3); hold off;
	semplot(t,reshape(mc(:,1,:),[],length(t)),color);
	semplot(t,reshape(mc(:,2:end,:),[],length(t)),'k');
	ylim([-0.012 0.02]);
	PlotHVLines([0 0.2],'v','k--');
	try
		xx = Restrict(sig{i,j},[-0.5 0.5]);
		plot(xx,repmat(max(ylim)-[0.05]*diff(ylim),size(xx))','color',color,'linewidth',5);
	catch
		drawnow
		[~,stats{i,j}] = CompareDistributionsNan(squeeze(mc(:,1,:)),reshape(squeeze(mc(:,2:end,:)),[],length(t)));
		try  sig{i,j} = t(FindInterval(stats{i,j}.above));
			xx = Restrict(sig{i,j},[-0.5 0.5]);
			plot(xx,repmat(max(ylim)-[0.05]*diff(ylim),size(xx))','color',color,'linewidth',5);
		end
	end
	if ~isempty(xx), titl = [titl ', sig from ' num2str(xx(1)) ' to ' num2str(xx(2))]; end
	xlabel('time from ripple (s)');
	title(titl);
end


%% Compute correlations between ripple activity and following corical activity (small window, deltaSpikes only)

movingWindow = -0.5:0.02:0.5;
t = movingWindow;
nSessions = size(X,1);
duration = 0.2;
nShifts = 5;

silenceWindow = [-0.015 0.015];
restrictToLocalDeltas = 1;
[mc2,mp2,mi2,mc3,mp3,mi3] = deal([]);
for j=3:4
	for ses=1:nSessions,
		[deltas,occ,groups,pfcUnit,ripples,pfc,hpc,~,session] = X{ses,:};
		[~,sessionName] = fileparts(session);
		occ = occ>0.5;
		groups = round(groups); if min(groups)>8,groups=groups-8;end
		pfcUnit =pfcUnit(:,1); if min(pfcUnit)>8,pfcUnit=pfcUnit-8;end
		if ~ismember(str2double(sessionName(4:6)),[199 211 214]),rank = tetrodeRank; else rank = tetrodeRank2;end
		pfcUnit = rank(pfcUnit); groups = rank(groups);
	
		[sessionPath,sessionName] = fileparts(session);
		occ = occ>0.5;
		groups = round(groups); if min(groups)>8,groups=groups-8;end
		peaks = X{ses,7};
		pfcUnit = pfcUnit(:,1); if min(pfcUnit)>8,pfcUnit=pfcUnit-8;end
		positions = dlmread(fullfile(fileparts(session),'positions'));
		behaviour = positions([1 end],1)';
		ripples = Restrict(ripples,[behaviour(end) Inf]);
		preceding = InIntervals(ripples(:,2),[deltas(:,2)-0.2 deltas(:,2)]);
		
		switch j
			case 3
				ok = ~preceding;
			case 4
				ok = preceding;
		end
		
		if sum(ok)<20 | sum(~ok)<20,
			nCorrs(ses,j) = 0;nEvents(ses,j) = 0;continue;
		end
		
		ripples = ripples(ok,:);
		rIntervals = [ripples(:,2)-duration/2 ripples(:,2)+duration/2]; % special: make them 200 ms for a change
		
		switch j
			case 3,
				source = CountSpikesIntervals(hpc,rIntervals);
				intervals = bsxfun(@plus,ripples(:,2)+0.13,silenceWindow);
				target = CountSpikesIntervals(pfc,intervals);
			case 4
				deltaSpikes = CountSpikesIntervals(pfc,bsxfun(@plus,deltas(:,2),silenceWindow));
				%     % Restrict to local deltas:
				if restrictToLocalDeltas,
					deltaSpikes(~occ(:,FindClosest(groups,pfcUnit))) = nan;
				end
				idx = FindClosest(deltas(:,2),ripples(:,2)+0.13);
				target = deltaSpikes(idx,:);
				intervals = bsxfun(@plus,deltas(FindClosest(deltas(:,2),ripples(:,2)+0.13),2),[-0.2 0]);
				source = CountSpikesIntervals(hpc,intervals);
				[~,good] = unique(idx);
				source = source(good,:); target = target(good,:);
		end
		cc = []; pp = [];
		thisc = []; thisp = [];
		for shift = 0:nShifts,
			[thiscShifted,thispShifted] = nancorr(target((1:end-nShifts)+shift,:),source((1:end-nShifts),:),'type','spearman'); thisc = [thisc thiscShifted(:)]; thisp = [thisp thispShifted(:)];
		end
		cc = cat(3,cc,thisc); pp = cat(3,pp,thisp);
		[idSource,idTarget] = meshgrid(1:size(source,2),1:size(target,2));
		switch j
			case 3,
				mc2 = cat(1,mc2,cc); mp2 = cat(1,mp2,pp); mi2 = cat(1,mi2,[idSource(:) idTarget(:)]);
				case 4,
				mc3 = cat(1,mc3,cc); mp3 = cat(1,mp3,pp); mi3 = cat(1,mi3,[idSource(:) idTarget(:)]);
		end
		nCorrs(ses,j) = numel(idSource);
		nEvents(ses,j) = size(source,1);
% 		ses
	end
end

%% Plot results

figure(2);
subplot(1,3,1);
g = Group(mc2(:,1),mc2(:,2:end));
anovabar(g(:,1),g(:,2));
set(gca,'xtick',1:2,'xticklabel',{'data (no delta)','shifted control'});
ylabel('mean corr between [-100 100]ms ripple window and [115 145]ms cortical spikes');
title(['No delta, p=' num2str(ranksum(g(g(:,2)==1),g(g(:,2)==2)))]);

subplot(1,3,2);
g = Group(mc3(:,1),mc3(:,2:end));
anovabar(g(:,1),g(:,2));
set(gca,'xtick',1:2,'xticklabel',{'data (delta)','shifted control'});
title(['Delta, p=' num2str(ranksum(g(g(:,2)==1),g(g(:,2)==2)))]);
ylabel('mean corr between [-100 100]ms ripple window and cortical deltaspikes');

subplot(1,3,3);
g = Group(mc2(:,1),mc3(:,1));
anovabar(g(:,1),g(:,2));
title(['No delta vs delta, p=' num2str(ranksum(g(g(:,2)==1),g(g(:,2)==2)))]);
set(gca,'xtick',1:2,'xticklabel',{'no delta','delta'});

EquateScales;


