function [correl]=multiple_correlation(data,tw,Time,znorm,shoulder,imagecolour,autonorm,ROI,PW50)

load(data.name);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

signals=whos('sig*');
if isempty(signals)==1;
    signals=(whos('*ch*'));
end
units = cell(1,length(signals));
for i =1:length(signals);
    cmd = char(strcat('units{i}=',signals(i).name));
    eval(cmd);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[y]=spike_selection(units,Time)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:length(units);

    for h=1:length(units);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    
        c=CORE_pxcorr(y{i},y{h},tw);
        correl{i}(h,:)=c;
               
        MeanCorrel{i}(h)=mean(c);
        standdev=std(c);
        confint{i}(h,1)=MeanCorrel{i}(h)-1.96*standdev/sqrt(length(c));
        confint{i}(h,2)=MeanCorrel{i}(h)+1.96*standdev/sqrt(length(c));
        clear c
        
    end    
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%Calculate peak%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[POffset,peakc, peakc1, HalfPeak, HPW,widths1,widths2,kindex]=peak_correl(correl,confint,signals,shoulder,MeanCorrel,PW50);
if znorm==1;
[zcorrel,peakn,confint,HalfPeak,HPW,widths1,widths2]=norm_corr(correl,shoulder,signals)
end
if autonorm==1;
[autopeakc]=auto_norm_corr(correl)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%Peak Offset%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
POffsetHist1={};POffsetHist2={};POffsetHist3={};
for i=1:(ROI-1);
    POffsetHist1{i}=POffset(i+1:ROI,i)';
end
for i=1:ROI;
    POffsetHist3{i}=POffset(ROI+1:end,i)';
end
for i=(ROI+1):length(POffset)-1;
    POffsetHist2{i-ROI}=POffset(i+1:end,i)';
end
POffsetHist1m=cell2mat(POffsetHist1);POffsetHist2m=cell2mat(POffsetHist2);POffsetHist3m=cell2mat(POffsetHist3);
binsize=1  %%%%%change this to alter the bin width of the peak offset histogram

figure(1);
set(gcf,'Name','Peak Offset Histogram ROI 1 (intra-correlation)');clf
hist(POffsetHist1m,-1000:binsize:1000)
ylabel('No. unit pairs/bin')
xlabel('Offset (ms)')
xlim([-PW50/2 PW50/2])
hdt=datacursormode
set(hdt,'UpdateFcn',{@pairinfo,signals,POffsetHist3,binsize,ROI,0})

figure(2);
set(gcf,'Name','Peak Offset Histogram ROI 2 (intra-correlation) ');clf
hist(POffsetHist2m,-1000:binsize:1000)
ylabel('No. unit pairs (/bin)')
xlabel('Offset (ms)')
xlim([-PW50/2 PW50/2])
hdt=datacursormode
set(hdt,'UpdateFcn',{@pairinfo,signals,POffsetHist3,binsize,ROI,1})

figure(3);
set(gcf,'Name','Peak Offset Histogram ROI 3 (inter-region correlation)');clf
hist(POffsetHist3m,-1000:binsize:1000)
ylabel('No. unit pairs/bin') 
xlabel('Offset (ms)')
xlim([-PW50/2 PW50/2])
hdt=datacursormode
set(hdt,'UpdateFcn',{@pairinfo,signals,POffsetHist3,binsize,ROI,1})
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%Bar plot for Half-Peak Width (HPW - PW50) 
R1={};R2={};R3={};
for i=1:ROI;
    R1{i}=HPW(find(isnan(HPW(i,i+1:ROI))==0));
    R2{i}=HPW(find(isnan(HPW(i+ROI,i+(ROI+1):end))==0));
    R3{i}=HPW(find(isnan(HPW(i,ROI+1:end))==0));
end
R1=cell2mat(R1);R2=cell2mat(R2);R3=cell2mat(R3);
Rmean(1)=mean(R1);Rmean(2)=mean(R2);Rmean(3)=mean(R3);
Rerror(1)=std(R1)/sqrt(length(R1));Rerror(2)=std(R2)/sqrt(length(R2));Rerror(3)=std(R3)/sqrt(length(R3));
figure(4);
set(gcf,'Name','Half-Peak Width (PW50)');clf
bar(Rmean,'k')
hold on
h=errorbar(Rmean,Rerror)
set(h(1),'color','k')
set(h,'LineStyle','none')
hold off
ylabel('Mean Width (PW50 ms)')
set(gca,'xticklabel',{['ROI 1'],['ROI 2'],['ROI 3']})

%%Bar Plot of kindex
R1={};R2={};R3={};
for i=1:ROI;
    R1{i}=kindex(i,i+1:ROI);
    R2{i}=kindex(i+ROI,i+(ROI+1):end);
    R3{i}=kindex(i,ROI+1:end);
end
R1=cell2mat(R1);R2=cell2mat(R2);R3=cell2mat(R3);
Rmean(1)=mean(R1);Rmean(2)=mean(R2);Rmean(3)=mean(R3);
Rerror(1)=std(R1)/sqrt(length(R1));Rerror(2)=std(R2)/sqrt(length(R2));Rerror(3)=std(R3)/sqrt(length(R3));
figure(5);
set(gcf,'Name','K-index');clf
bar(Rmean,'k')
hold on
h=errorbar(Rmean,Rerror)
set(h(1),'color','k')
set(h,'LineStyle','none')
hold off
ylabel('Mean K-index')
set(gca,'xticklabel',{['ROI 1'],['ROI 2'],['ROI 3']})

%%%Plot Grid%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if exist('zcorrel')==0 & autonorm==1
plot_grid(signals,imagecolour,autonorm,autopeakc)
hdt = datacursormode;
set(hdt,'UpdateFcn',{@correl_hist,correl,tw,znorm,signals,MeanCorrel,confint, HalfPeak, HPW,widths1,widths2,kindex})
datacursormode on
elseif exist('zcorrel')==0 & autonorm==2
plot_grid(signals,imagecolour,autonorm,peakc,peakc1)
hdt = datacursormode;
set(hdt,'UpdateFcn',{@correl_hist,correl,tw,znorm,signals,MeanCorrel,confint, HalfPeak, HPW,widths1,widths2,kindex})
%set(hdt,'windowbuttonmotionfcn',@mousemove)
else
plot_grid(signals,imagecolour,autonorm,peakn)
hdt = datacursormode;
set(hdt,'UpdateFcn',{@correl_hist,zcorrel,tw,znorm,signals,MeanCorrel,confint,HalfPeak,HPW,widths1,widths2,kindex})
end








