function [] = figure6;

%
% function [] = figure6;
%
% plots figure 6 of 
%
% Einhäuser, W., Thomassen, S., & Bendixen, A. (2017). Using binocular
% rivalry to tag foreground sounds: towards an objective visual
% measure for auditory multistability. Journal of Vision, 17:34,
% 1-19. 
%
%

load exp1_data gain ixLow ixHigh ixBoth

[meanTraceHighToLow(1:8,:),...
	  meanTraceLowToHigh(1:8,:),...
	  meanTraceHighToBoth(1:8,:),...
	  meanTraceBothToHigh(1:8,:),...
	  meanTraceLowToBoth(1:8,:),...
	  meanTraceBothToLow(1:8,:)] = computeAlignedTraces(ixLow,ixHigh,ixBoth,gain,[5:8,10:13]);


load exp2_data gain ixLow ixHigh ixBoth

[meanTraceHighToLow(9:16,:),...
	  meanTraceLowToHigh(9:16,:),...
	  meanTraceHighToBoth(9:16,:),...
	  meanTraceBothToHigh(9:16,:),...
	  meanTraceLowToBoth(9:16,:),...
	  meanTraceBothToLow(9:16,:)] = ...
	computeAlignedTraces(ixLow,ixHigh,ixBoth,gain,[10:13,15:18]);



load exp3_data gain ixLow ixHigh ixBoth

[meanTraceHighToLow(17:24,:),...
	  meanTraceLowToHigh(17:24,:),...
	  meanTraceHighToBoth(17:24,:),...
	  meanTraceBothToHigh(17:24,:),...
	  meanTraceLowToBoth(17:24,:),...
	  meanTraceBothToLow(17:24,:)] = computeAlignedTraces(ixLow,ixHigh,ixBoth,gain,[5:8,10:13]);


figure(6);
clf;

subplot(2,3,1);
plot([-1999:2000]/500,squeeze(nanmean(meanTraceHighToLow)),'b');
hold on
plot([-1999:2000]/500,squeeze(nanmean(meanTraceLowToHigh)),'r');
plot([-1999:2000]/500,squeeze(nanmean(meanTraceLowToHigh)+nanstd(meanTraceLowToHigh)./sqrt(sum(isfinite(meanTraceLowToHigh)))),'m');
plot([-1999:2000]/500,squeeze(nanmean(meanTraceLowToHigh)-nanstd(meanTraceLowToHigh)./sqrt(sum(isfinite(meanTraceLowToHigh)))),'m');
plot([-1999:2000]/500,squeeze(nanmean(meanTraceHighToLow)+nanstd(meanTraceHighToLow)./sqrt(sum(isfinite(meanTraceHighToLow)))),'c');
plot([-1999:2000]/500,squeeze(nanmean(meanTraceHighToLow)-nanstd(meanTraceHighToLow)./sqrt(sum(isfinite(meanTraceHighToLow)))),'c');
ylim([-.4 .4]);
ylabel('gain');
D = meanTraceLowToHigh-meanTraceHighToLow;
[h,p] = ttest(meanTraceLowToHigh,meanTraceHighToLow);
subplot(2,3,4);
plot([-1999:2000]/500,-log10(p));
xlabel('time to report [s]');
ylabel('-log_{10}(p)');
subplot(2,3,2);
plot([-1999:2000]/500,squeeze(nanmean(meanTraceHighToBoth)),'b');
hold on
plot([-1999:2000]/500,squeeze(nanmean(meanTraceLowToBoth)),'r');
plot([-1999:2000]/500,squeeze(nanmean(meanTraceLowToBoth)+nanstd(meanTraceLowToBoth)./sqrt(sum(isfinite(meanTraceLowToBoth)))),'m');
plot([-1999:2000]/500,squeeze(nanmean(meanTraceLowToBoth)-nanstd(meanTraceLowToBoth)./sqrt(sum(isfinite(meanTraceLowToBoth)))),'m');
plot([-1999:2000]/500,squeeze(nanmean(meanTraceHighToBoth)+nanstd(meanTraceHighToBoth)./sqrt(sum(isfinite(meanTraceHighToBoth)))),'c');
plot([-1999:2000]/500,squeeze(nanmean(meanTraceHighToBoth)-nanstd(meanTraceHighToBoth)./sqrt(sum(isfinite(meanTraceHighToBoth)))),'c');
ylim([-.4 .4]);
ylabel('gain');
D = meanTraceLowToBoth-meanTraceHighToBoth;
[h,p] = ttest(meanTraceLowToBoth,meanTraceHighToBoth);
subplot(2,3,5);
plot([-1999:2000]/500,-log10(p));
xlabel('time to report [s]');
ylabel('-log_{10}(p)');

subplot(2,3,3);
plot([-1999:2000]/500,squeeze(nanmean(meanTraceBothToLow)),'b');
hold on
plot([-1999:2000]/500,squeeze(nanmean(meanTraceBothToHigh)),'r');
plot([-1999:2000]/500,squeeze(nanmean(meanTraceBothToHigh)+nanstd(meanTraceBothToHigh)./sqrt(sum(isfinite(meanTraceBothToHigh)))),'m');
plot([-1999:2000]/500,squeeze(nanmean(meanTraceBothToHigh)-nanstd(meanTraceBothToHigh)./sqrt(sum(isfinite(meanTraceBothToHigh)))),'m');
plot([-1999:2000]/500,squeeze(nanmean(meanTraceBothToLow)+nanstd(meanTraceBothToLow)./sqrt(sum(isfinite(meanTraceBothToLow)))),'c');
plot([-1999:2000]/500,squeeze(nanmean(meanTraceBothToLow)-nanstd(meanTraceBothToLow)./sqrt(sum(isfinite(meanTraceBothToLow)))),'c');
ylim([-.4 .4]);
ylabel('gain');
D = meanTraceBothToHigh-meanTraceBothToLow;
[h,p] = ttest(meanTraceBothToHigh,meanTraceBothToLow);
subplot(2,3,6);
plot([-1999:2000]/500,-log10(p));
xlabel('time to report [s]');
ylabel('-log_{10}(p)');








function [meanTraceHighToLow,...
	  meanTraceLowToHigh,...
	  meanTraceHighToBoth,...
	  meanTraceBothToHigh,...
	  meanTraceLowToBoth,...
	  meanTraceBothToLow] = computeAlignedTraces(idxLowTone,idxHighTone,idxBoth,gain,blocks)


nSamples = 2000;

meanTraceHighToLow = nan(8,2*nSamples);
meanTraceLowToHigh = nan(8,2*nSamples);

meanTraceHighToBoth = nan(8,2*nSamples);
meanTraceBothToHigh = nan(8,2*nSamples);

meanTraceBothToLow = nan(8,2*nSamples);
meanTraceLowToBoth = nan(8,2*nSamples);

transCount=zeros(1,8);

for subNum = 1:8

  trace_lowToHigh   = nan(0,2*nSamples);
  trace_highToLow  =  nan(0,2*nSamples);
  
  trace_bothToHigh   = nan(0,2*nSamples);
  trace_highToBoth  =  nan(0,2*nSamples);
  
  trace_bothToLow   = nan(0,2*nSamples);
  trace_lowToBoth  =  nan(0,2*nSamples);
  
  for block = blocks
	% extract onsets and offsets
	
	if ~isempty(idxLowTone{subNum,block})
	  ixx           = find(diff(idxLowTone{subNum,block})>1);
	  onsetLowTone  = idxLowTone{subNum,block}([1,ixx+1]);
	  offsetLowTone = idxLowTone{subNum,block}([ixx,length(idxLowTone)]);
	end
	
	if ~isempty(idxHighTone{subNum,block})
	  ixx           = find(diff(idxHighTone{subNum,block})>1);
	  onsetHighTone  = idxHighTone{subNum,block}([1,ixx+1]);
	  offsetHighTone = idxHighTone{subNum,block}([ixx,length(idxHighTone)]);
	end
	
	if ~isempty(idxBoth{subNum,block})
	  ixx           = find(diff(idxBoth{subNum,block})>1);
	  onsetInt  = idxBoth{subNum,block}([1,ixx+1]);
	  offsetInt = idxBoth{subNum,block}([ixx,length(idxBoth)]);
	end
	
	

	% now search high tone onsets that are within 100ms of a low tone
	% offset and vice versa
	for k = 1:length(onsetHighTone)
	  currentT = onsetHighTone(k);
	  ixLowOff = find((currentT>offsetLowTone)&((currentT-offsetLowTone)<50));
  
	  if ~isempty(ixLowOff)
		lowToHigh=currentT;
		lowToHigh_start=onsetLowTone(ixLowOff); % onset of the low tone
		lowToHigh_end  =offsetHighTone(k);
		
		Lpos = min(nSamples,lowToHigh_end-lowToHigh);
		Lneg = min(nSamples,lowToHigh-lowToHigh_start);
	
		tmp = nan(1,2*nSamples);
		tmp((nSamples-Lneg+1):nSamples)  = gain{subNum,block,1}((lowToHigh-Lneg+1):lowToHigh);
		tmp((nSamples+1):(nSamples+Lpos))  = gain{subNum,block,1}((lowToHigh+1):(lowToHigh+Lpos));
		trace_lowToHigh(end+1,:) = tmp;
		
	  end
	  
	end
	
	for k = 1:length(onsetLowTone)
	  currentT = onsetLowTone(k);
	  ixHighOff = find((currentT>offsetHighTone)&((currentT-offsetHighTone)<50));
	  
	  if ~isempty(ixHighOff)
		highToLow=currentT;
		highToLow_start=onsetHighTone(ixHighOff); % onset of the high tone
		highToLow_end  =offsetLowTone(k);
		
		Lpos = min(nSamples,highToLow_end-highToLow);
		Lneg = min(nSamples,highToLow-highToLow_start);
		
		tmp = nan(1,2*nSamples);
		tmp((nSamples-Lneg+1):nSamples)  = gain{subNum,block,1}((highToLow-Lneg+1):highToLow);
		tmp((nSamples+1):(nSamples+Lpos))  = gain{subNum,block,1}((highToLow+1):(highToLow+Lpos));
		trace_highToLow(end+1,:) = tmp;
	  end
	  
	end
	
	
	for k = 1:length(onsetInt)
	  currentT = onsetInt(k);
	  ixLowOff = find((currentT>offsetLowTone)&((currentT-offsetLowTone)<50));
	  
	  if ~isempty(ixLowOff)
		lowToBoth=currentT;
		lowToBoth_start=onsetLowTone(ixLowOff); % onset of the low tone
		lowToBoth_end  =offsetInt(k);
		
		Lpos = min(nSamples,lowToBoth_end-lowToBoth);
		Lneg = min(nSamples,lowToBoth-lowToBoth_start);
		
		tmp = nan(1,2*nSamples);
		tmp((nSamples-Lneg+1):nSamples)  = gain{subNum,block,1}((lowToBoth-Lneg+1):lowToBoth);
		tmp((nSamples+1):(nSamples+Lpos))  = gain{subNum,block,1}((lowToBoth+1):(lowToBoth+Lpos));
		trace_lowToBoth(end+1,:) = tmp;
		
	  end
  
	end
	
	for k = 1:length(onsetInt)
	  currentT = onsetInt(k);
	  ixHighOff = find((currentT>offsetHighTone)&((currentT-offsetHighTone)<50));
	  
	  if ~isempty(ixHighOff)
		highToBoth=currentT;
		highToBoth_start=onsetHighTone(ixHighOff); % onset of the high tone
		highToBoth_end  =offsetInt(k);
		
		Lpos = min(nSamples,highToBoth_end-highToBoth);
		Lneg = min(nSamples,highToBoth-highToBoth_start);
		
		tmp = nan(1,2*nSamples);
		tmp((nSamples-Lneg+1):nSamples)  = gain{subNum,block,1}((highToBoth-Lneg+1):highToBoth);
		tmp((nSamples+1):(nSamples+Lpos))  = gain{subNum,block,1}((highToBoth+1):(highToBoth+Lpos));
		trace_highToBoth(end+1,:) = tmp;
		
	  end
  
	end
	
	
	
	for k = 1:length(onsetLowTone)
	  currentT = onsetLowTone(k);
	  ixBothOff = find((currentT>offsetInt)&((currentT-offsetInt)<50)&((currentT-onsetInt)>50));
	  
	  if ~isempty(ixBothOff)
		bothToLow=currentT;
		bothToLow_start=onsetInt(ixBothOff); % onset of the high tone
		bothToLow_end  =offsetLowTone(k);
		
		Lpos = min(nSamples,bothToLow_end-bothToLow);
		Lneg = min(nSamples,bothToLow-bothToLow_start);
		
		tmp = nan(1,2*nSamples);
		tmp((nSamples-Lneg+1):nSamples)  = gain{subNum,block,1}((bothToLow-Lneg+1):bothToLow);
		tmp((nSamples+1):(nSamples+Lpos))  = gain{subNum,block,1}((bothToLow+1):(bothToLow+Lpos));
		trace_bothToLow(end+1,:) = tmp;
	  end
	  
	end
	
	for k = 1:length(onsetHighTone)
	  currentT = onsetHighTone(k);
	  ixBothOff = find((currentT>offsetInt)&((currentT-offsetInt)<50)&((currentT-onsetInt)>50));
	  
	  if ~isempty(ixBothOff)
		bothToHigh=currentT;
		bothToHigh_start=onsetInt(ixBothOff); % onset of the high tone
		bothToHigh_end  =offsetHighTone(k);
		
		Lpos = min(nSamples,bothToHigh_end-bothToHigh);
		Lneg = min(nSamples,bothToHigh-bothToHigh_start);
		
		tmp = nan(1,2*nSamples);
		tmp((nSamples-Lneg+1):nSamples)  = gain{subNum,block,1}((bothToHigh-Lneg+1):bothToHigh);
		tmp((nSamples+1):(nSamples+Lpos))  = gain{subNum,block,1}((bothToHigh+1):(bothToHigh+Lpos));
		trace_bothToHigh(end+1,:) = tmp;
	  end
  
	end

  end	


  meanTraceHighToLow(subNum,:) = nanmean(trace_highToLow);
  meanTraceLowToHigh(subNum,:) = nanmean(trace_lowToHigh);
  
  meanTraceHighToBoth(subNum,:) = nanmean(trace_highToBoth);
  meanTraceBothToHigh(subNum,:) = nanmean(trace_bothToHigh);

  meanTraceLowToBoth(subNum,:) = nanmean(trace_lowToBoth);
  meanTraceBothToLow(subNum,:) = nanmean(trace_bothToLow);

    
end


