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

%% Compute correlations

movingWindow = -1:0.01:1;
silenceWindow = [-0.015 0.015];
t = movingWindow;
nSessions = size(X,1);
[mc,mp,mi,gains,gainsMUA,shGains,shGainsMUA] = deal([]); % m is for moving
[nSpikes,nCorr] = deal(nan(nSessions,1));
duration = 0.2;

tetrodeAnatomy = [4 5 3 6 2 7 1 8]; % from anterior to posterior
tetrodeAnatomy2 = [8 4 7 3 6 2 5 1]; % from anterior to posterior; for rats 211, 199, and 214
[~,tetrodeRank] = sort(tetrodeAnatomy); % for each tetrode, its rank in the anatomical anterio-posterior axis
[~,tetrodeRank2] = sort(tetrodeAnatomy2); % for each tetrode, its rank in the anatomical anterio-posterior axis; for rats 211, 199, and 214
meanDeltaRate = [];
nShifts = 3;
for ses=1:nSessions,
	[deltas,occ,groups,pfcUnit,ripples,pfc,hpc,varargin,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
	
	prediction = nan(size(occ));
	rIntervals = [ripples(:,2)-0.05 ripples(:,2)+0.05];
	if ~ismember(str2double(sessionName(4:6)),[199 211 214]),rank = tetrodeRank; else rank = tetrodeRank2;end
	pfcUnit = rank(pfcUnit); groups = rank(groups);
	source = CountSpikesIntervals(pfc,bsxfun(@plus,deltas(:,2),silenceWindow));
	meanDeltaRate = [meanDeltaRate; nanmean(source)'];
	
	cc = []; pp = [];
	for i=1:length(movingWindow),
		intervals = bsxfun(@plus,deltas(:,2),[-1 1]*duration/2+movingWindow(i));
		target = CountSpikesIntervals(hpc,intervals);
		% restrict to ripples,
		ripplesIn = IntervalsIntersect(intervals,rIntervals);
		target(~ripplesIn,:) = nan;
		thisc = []; thisp = [];
		for shift = 0:nShifts,
			[thiscShifted,thispShifted] = nancorr(source((1:end-nShifts)+shift,:),target((1:end-nShifts),:),'type','spearman'); thisc = [thisc thiscShifted(:)]; thisp = [thisp thispShifted(:)];
		end
		cc = cat(3,cc,thisc); pp = cat(3,pp,thisp);
	end
	[idTarget,idSource] = meshgrid(1:size(target,2),1:size(source,2));
	mc = cat(1,mc,cc); mp = cat(1,mp,pp); mi = cat(1,mi,[idTarget(:) idSource(:)]);
	nCorr(ses,1) = numel(idTarget);
	ses
end

sesh = repelem((1:size(X,1))',nCorr); clear day sync
for i=1:size(X,1),
	if strcmp(X{i,16}{8},'Test'),day(i,1) = 2; else; day(i,1) = 1; end
	if strcmp(X{i,16}{6},'CTRL'),sync(i,1) = 0; else; sync(i,1) = 1; end
end

t = movingWindow;
%%

figure(1)
smooth=1;
subplot(2,1,1); cla
hold off
semplot(t,squeeze(reshape(mc(:,2:end,:),[],size(mc,3))),'k',smooth);
semplot(movingWindow,squeeze(mc(:,1,:)),[0.7 0.3 0],smooth);
PlotIntervals(silenceWindow,'color',[0.5 0.5 0.5],'ylim',min(ylim)+[0 0.05]*diff(ylim));
ylim([-0.008 0.015]);

PlotHVLines(0,'v','k--');
ylabel('average correlation');
xlabel(['center of the ' num2str(diff(silenceWindow)*1000) 'ms delta window']);
title(['2B for ' num2str(diff(silenceWindow)*1000)  'ms window around delta wave']);

drawnow
[~,stats]  = CompareDistributionsNan(squeeze(mc(:,1,:)),reshape(squeeze(mc(:,2:end,:)),[],length(t)));
PlotIntervals(t(FindInterval(stats.above)),'color',[1 0.5 0.5]);
PlotIntervals(t(FindInterval(stats.below)),'color',[0.5 0.5 1]);
xlim([-1 1]*0.5);

subplot(2,1,2);
smooth = [4];
[h1,ht] = hist(squeeze(mc(:,1,:)),linspace(-0.3,0.3,50));
[h0,ht] = hist((squeeze(mc(:,2,:))),ht); %h0 = repmat(h0(:,1),1,size(h1,2))./size(h1,2);

PlotColorMap(100*Smooth(bsxfun(@rdivide,(h1-h0),sum(h1+h0)./2),smooth,'type','cc'),'x',t,'y',ht);
PlotHVLines(0,'v','k--'); PlotHVLines(0,'h','k--');
xlim([-1 1]*0.5);
clim([-1 1]*0.2); clabel('Enrichment (%)');
axis('square');







