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

%% 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));
	
	% Restrict to local deltas:
	source(~occ(:,FindClosest(groups,pfcUnit))) = nan;
	
	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(:)];
			
			% GLM
			if shift==0 && i==FindClosest(movingWindow,-0.1);
				nUnits = size(source,2); units(ses,1) = nUnits;
				[gain,gainMUA] = deal(nan(nUnits,1));
				[shGain,shGainMUA] = deal(nan(nUnits,1000));
				
				for unit = 1:nUnits,
					s = target; t = source(:,unit);
					nans = isnan(t) | isnan(s(:,1)); % events that were not detected as delta waves on this channel
					s(nans,:) = []; t(nans) = [];
					if sum(t>0)>=5,
						[gain(unit,1),shGain(unit,:),er] = GLMgain(s,double(t>0),'nsets',sum(t>0),'dist','normal','link','identity','minevents',4);
						[gainMUA(unit,1),shGainMUA(unit,:)] = GLMgain(sum(s,2),double(t>0),'nsets',sum(t>0),'dist','normal','link','identity','minevents',5);
					end
				end
				nSpikes = [nSpikes; nansum(target>0)'];
				gains = [gains;gain];
				shGains = [shGains;shGain];
				gainsMUA = [gainsMUA;gainMUA];
				shGainsMUA = [shGainsMUA;shGainMUA];
			end
		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 rat
ratNumber = @(x) str2double(x((3:5)+min(strfind(lower(x),'rat'))));
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
	session = X{i,10};
	rat(i,1) = ratNumber(session);
end

t = movingWindow;
%%

ratID = Renumber(rat(sesh)); % rename rats as 1,2,3...
uniqueRat = rat(sesh);
figure(2)
smooth=1;
for i=1:15,
	subplot(5,3,i); cla
	ok = sesh~=i;
	hold off
	semplot(t,squeeze(reshape(mc(ok,2:end,:),[],size(mc,3))),'k',smooth);
	semplot(movingWindow,squeeze(mc(ok,1,:)),[0.7 0.3 0],smooth);
	ylim([-0.008 0.015]);
	
	PlotIntervals(silenceWindow,'color',[0.5 0.5 0.5]/2,'ylim',min(ylim)+[0 0.1]*diff(ylim));
	PlotHVLines(0,'v','k--');PlotHVLines(0,'h','k--');
	if rem(i,3)==1, % i.e. first column
		ylabel('average correlation');
	end
	if i>12, % i.e. last row
		xlabel(['center of the ' num2str(diff(silenceWindow)*1000) 'ms delta window']);
	end
	
	try
		plot(t(FindInterval(sesStats{i,1}.above))',repmat(max(ylim)-[0.05]*diff(ylim),size(FindInterval(sesStats{i,1}.above)))','color',[0.7 0.3 0],'linewidth',5);
	catch
		drawnow
		[~,sesStats{i,1}]  = CompareDistributionsNan(squeeze(mc(:,1,:)),reshape(squeeze(mc(:,2:end,:)),[],length(t)));
		plot(t(FindInterval(sesStats{i,1}.above))',repmat(max(ylim)-[0.05]*diff(ylim),size(FindInterval(sesStats{i,1}.above)))','color',[0.7 0.3 0],'linewidth',5);
	end
	xlim([-1 1]*0.5);
	title(['Without ' out2(@fileparts,X{i,10})]);
end

%%


ratID = Renumber(rat(sesh)); % rename rats as 1,2,3...
uniqueRat = rat(sesh);
figure(2)
smooth=1;
for i=1:9,
	subplot(3,3,i); %cla
	ok = ratID~=i;
	hold off
	semplot(t,squeeze(reshape(mc(ok,2:end,:),[],size(mc,3))),'k',smooth);
	semplot(movingWindow,squeeze(mc(ok,1,:)),[0.7 0.3 0],smooth);
	ylim([-0.008 0.015]);
	
	PlotIntervals(silenceWindow,'color',[0.5 0.5 0.5]/2,'ylim',min(ylim)+[0 0.1]*diff(ylim));
	PlotHVLines(0,'v','k--');PlotHVLines(0,'h','k--');
	if rem(i,3)==1, % i.e. first column
		ylabel('average correlation');
	end
	if i>6, % i.e. last row
		xlabel(['center of the ' num2str(diff(silenceWindow)*1000) 'ms delta window']);
	end
	drawnow
	[~,ratStats{i,1}]  = CompareDistributionsNan(squeeze(mc(:,1,:)),reshape(squeeze(mc(:,2:end,:)),[],length(t)));
	plot(t(FindInterval(ratStats{i,1}.above))',repmat(max(ylim)-[0.05]*diff(ylim),size(FindInterval(ratStats{i,1}.above)))','color',[0.7 0.3 0],'linewidth',5);
	xlim([-1 1]*0.5);
	title(['Leaving out rat ' num2str(rat(sesh(find(ratID==i,1)))) ]);
	if sum(~ok)==0,
		title('All rats');
	end
end
