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,~,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);
	
	deltaSpikes = CountSpikesIntervals(pfc,bsxfun(@plus,deltas(:,2),silenceWindow));
	meanDeltaRate = [meanDeltaRate; nanmean(deltaSpikes)'];
	
	cc = []; pp = [];
	for i=1:length(movingWindow),
		intervals = bsxfun(@plus,deltas(:,2),[-1 1]*duration/2+movingWindow(i));
		target = CountSpikesIntervals(hpc,intervals);
		ripplesIn = IntervalsIntersect(intervals,rIntervals); % restrict to ripples
		target(~ripplesIn,:) = nan;
		
		thisc = []; thisp = [];
		for shift = 0:nShifts,
			[thiscShifted,thispShifted] = nancorr(deltaSpikes((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(deltaSpikes,2); units(ses,1) = nUnits;
				[gain,gainMUA] = deal(nan(nUnits,1));
				[shGain,shGainMUA] = deal(nan(nUnits,1000));
				
				for unit = 1:nUnits,
					s = target; t = deltaSpikes(:,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(deltaSpikes,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;
mc00 = mc;

%% Ignore all non-sig correlations

mc(mp>0.05) = nan;
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(['Ignoring all non-significant correlations']);

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

subplot(2,1,2);
smooth = [4];
[h1,ht] = hist(squeeze(mc(:,1,51:151)),linspace(-0.3,0.3,50));
[h0,ht] = hist((squeeze(mc(:,3,51:151))),ht); %h0 = repmat(h0(:,1),1,size(h1,2))./size(h1,2);
% subplot(1,1,1);
PlotColorMap(100*Smooth(bsxfun(@rdivide,(h1-h0),sum(h1+h0)./2),smooth,'type','cc'),'x',t(51:151),'y',ht);
PlotHVLines(0,'v','k--'); PlotHVLines(0,'h','k--');
xlim([-1 1]*0.5);
clim([-1 1]*1.5);
clabel('Enrichment (%)');
axis('square');

%% Take only partner correlations (i.e. correlated within the right window)

[~,~,u(:,1)] = unique([sesh mi(:,1)],'rows');
[~,~,u(:,2)] = unique([sesh mi(:,2)],'rows');

sig = mp(:,1,FindClosest(t,-0.1))<0.05;
partnerCells = find((Accumulate(u(:,2),sig)>0));
good = ismember(u(:,2),partnerCells);

% Make a table for how many HPC partners each cell has
table = uniqueID; % session, number
table(:,3) = Accumulate(u(:,2),sig); % number of partners (between 0 and 4)
table(:,4) = Accumulate(u(:,2),~sig); % number of non-partners (min. 3 for partnered cells)

mc = mc00;
mc(~good,:,:) = nan;
figure(1)
smooth=1;
subplot(2,2,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(['Only for PFC cells with partners, n=' num2str(numel(partnerCells))]);

try
	plot(t(FindInterval(statsPartnerOnly.above))',repmat(max(ylim)-[0.05]*diff(ylim),size(FindInterval(statsPartnerOnly.above)))','color',[0.7 0.3 0],'linewidth',5);
catch
	drawnow
	[~,statsPartnerOnly]  = CompareDistributionsNan(squeeze(mc(:,1,:)),reshape(squeeze(mc(:,3:end,:)),[],length(t)));
	plot(t(FindInterval(statsPartnerOnly.above))',repmat(max(ylim)-[0.05]*diff(ylim),size(FindInterval(statsPartnerOnly.above)))','color',[0.7 0.3 0],'linewidth',5);
end
xlim([-1 1]*0.5);
axis('square');

subplot(2,2,3);
smooth = [4];
[h1,ht] = hist(squeeze(mc(:,1,51:151)),linspace(-0.3,0.3,50));
[h0,ht] = hist((squeeze(mc(:,3,51:151))),ht); %h0 = repmat(h0(:,1),1,size(h1,2))./size(h1,2);
% subplot(1,1,1);
PlotColorMap(100*Smooth(bsxfun(@rdivide,(h1-h0),sum(h1+h0)./2),smooth,'type','cc'),'x',t(51:151),'y',ht);
PlotHVLines(0,'v','k--'); PlotHVLines(0,'h','k--');
xlim([-1 1]*0.5);
clim([-1 1]*0.5);
clabel('Enrichment (%)');
axis('square');

subplot(2,2,2);
anovabar(([gainsMUA(partnerCells) gains(partnerCells)]-1)*100);
set(gca,'xticklabels',{['summed,p=' num2str(out2(@ttest2,gainsMUA(partnerCells),reshape(shGainsMUA(partnerCells),[],1)))], ['single units,p=' num2str(out2(@ttest2,gains(partnerCells),reshape(shGains(partnerCells,:),[],1)))]});
ylabel('Prediction gain (%)');
title(['p=' num2str(signrank(gainsMUA(partnerCells),gains(partnerCells))) ', between-group signrank']);


subplot(2,2,4);
% Try GLM-predicting partner cells via HPC partners or non-partners!
if ~exist('gainsYes','var');
[gainsYes,gainsNo,gainsNoMUA,shYes,shNo,shNoMUA] = deal([]);
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
	
	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);
	
	deltaSpikes = CountSpikesIntervals(pfc,bsxfun(@plus,deltas(:,2),silenceWindow));

	meanDeltaRate = [meanDeltaRate; nanmean(deltaSpikes)'];
	
	i=FindClosest(movingWindow,-0.1);
	
	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;
	
	nUnits = size(deltaSpikes,2); units(ses,1) = nUnits;
	[gain,gainMUA] = deal(nan(nUnits,1));
	[shGain,shGainMUA] = deal(nan(nUnits,1000));
	
	pValue = mp(sesh==ses,1,i);
	ids = mi(sesh==ses,:);
	partnered = unique(ids(pValue<0.05,2));
	[gainYes,gainNo,gainNoMUA] = deal(nan(nUnits,1));
	[sYes,sNo,sNoMUA] = deal(nan(nUnits,1000));
	for unit = partnered',
		s = target; t = deltaSpikes(:,unit);
		nans = isnan(t) | isnan(s(:,1)); % events that were not detected as delta waves on this channel
		s(nans,:) = []; t(nans) = [];
		partners = ids(pValue<0.05 & ids(:,2)==unit);
		if sum(t>0)>=5,
			[gainYes(unit,1),sYes(unit,:)] = GLMgain(s(:,partners),double(t>0),'nsets',sum(t>0),'dist','normal','link','identity','minevents',4);
			[gainNo(unit,1),sNo(unit,:)] = GLMgain(s(:,~ismember(1:size(s,2),partners)),double(t>0),'nsets',sum(t>0),'dist','normal','link','identity','minevents',5);
			[gainNoMUA(unit,1),sNoMUA(unit,:)] = GLMgain(sum(s(:,~ismember(1:size(s,2),partners)),2),double(t>0),'nsets',sum(t>0),'dist','normal','link','identity','minevents',5);
		end
	end
	nSpikes = [nSpikes; nansum(target>0)'];
	gainsYes = [gainsYes;gainYes];shYes = [shYes;sYes];
	gainsNo = [gainsNo;gainNo];shNo = [shNo;sNo];
	gainsNoMUA = [gainsNoMUA;gainNoMUA];shNoMUA = [shNoMUA;sNoMUA];
	
	ses
end
t = movingWindow;
end

anovabar(([gainsNoMUA gainsYes]-1)*100);
set(gca,'xticklabels',{['summed nonpartner,p=' num2str(out2(@ttest2,gainsNoMUA,shNoMUA(:)))], ['single partner units,p=' num2str(out2(@ttest2,gainsYes,shYes(:)))]});
title(['p=' num2str(signrank(gainsYes-gainsNoMUA)) ', between-group signrank']);
ylabel('Prediction gain (%)');

%% Take only PFC cells without partners (i.e. correlated within the right window)

[~,~,u(:,1)] = unique([sesh mi(:,1)],'rows');
[uniqueID,~,u(:,2)] = unique([sesh mi(:,2)],'rows');

sig = mp(:,1,FindClosest(t,-0.1))<0.05;
nonpartnerCells = find(~(Accumulate(u(:,2),sig)>0));
good = ismember(u(:,2),nonpartnerCells);

mc = mc00;
mc(~good,:,:) = nan;
figure(1)
smooth=1;
subplot(2,2,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(['Only for PFC cells without partners, n=' num2str(numel(nonpartnerCells))]);

try
	plot(t(FindInterval(statsNonPartnerOnly.above))',repmat(max(ylim)-[0.05]*diff(ylim),size(FindInterval(statsNonPartnerOnly.above)))','color',[0.7 0.3 0],'linewidth',5);
catch
	drawnow
	[~,statsNonPartnerOnly]  = CompareDistributionsNan(squeeze(mc(:,1,:)),reshape(squeeze(mc(:,3:end,:)),[],length(t)));
	plot(t(FindInterval(statsNonPartnerOnly.above))',repmat(max(ylim)-[0.05]*diff(ylim),size(FindInterval(statsNonPartnerOnly.above)))','color',[0.7 0.3 0],'linewidth',5);
end

xlim([-1 1]*0.5);
axis('square');

subplot(2,2,3);
smooth = [4];
[h1,ht] = hist(squeeze(mc(:,1,51:151)),linspace(-0.3,0.3,50));
[h0,ht] = hist((squeeze(mc(:,3,51:151))),ht); %h0 = repmat(h0(:,1),1,size(h1,2))./size(h1,2);
% subplot(1,1,1);
PlotColorMap(100*Smooth(bsxfun(@rdivide,(h1-h0),sum(h1+h0)./2),smooth,'type','cc'),'x',t(51:151),'y',ht);
PlotHVLines(0,'v','k--'); PlotHVLines(0,'h','k--');
xlim([-1 1]*0.5);
clim([-1 1]*0.5);
clabel('Enrichment (%)');
axis('square');

subplot(1,2,2);
anovabar(([gainsMUA(nonpartnerCells) gains(nonpartnerCells)]-1)*100);
set(gca,'xticklabels',{['summed,p=' num2str(out2(@ttest2,gainsMUA(nonpartnerCells),reshape(shGainsMUA(nonpartnerCells,:),[],1)))], ['single units,p=' num2str(out2(@ttest2,gains(nonpartnerCells),reshape(shGains(nonpartnerCells,:),[],1)))]});
title(['p=' num2str(ranksum(gainsMUA(nonpartnerCells),gains(nonpartnerCells))) ', between-group ranksum']);
ylabel('Prediction gain (%)');







