function [CS1,IS1,CS2,IS2,CSMat,ISMat,IndVec1,IndVec2]=MutualInfo(FuncMat,Base,NShotPerCycle,MinLag,MaxLag)

%	MutualInfo.m
%
% 	This function is designed to calculate the divergences 
%	between two statistical distributions: 
%	integrated square error (IS) and Cauchy-Schwartz measure  (CS) 

%   This software is protected by german copyright and international treaties.             
%   Copyright 2004 Qun Zhao, Markus Junghöfer & Peter Peyk. All Rights Reserved.                     
%                                                                                          
%   THIS SOFTWARE AND DOCUMENTATION IS PROVIDED "AS IS," AND COPYRIGHT HOLDERS MAKE        
%   NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO,    
%   WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR THAT THE USE    
%   OF THE SOFTWARE OR DOCUMENTATION WILL NOT INFRINGE ANY THIRD PARTY PATENTS, COPYRIGHTS,
%   TRADEMARKS OR OTHER RIGHTS. COPYRIGHT HOLDERS WILL NOT BE LIABLE FOR ANY DIRECT,       
%   INDIRECT, SPECIAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF ANY USE OF THE SOFTWARE      
%   OR DOCUMENTATION. 

if nargin<5; MaxLag=0; end
if nargin<4; MinLag=0; end
if nargin<3; NShotPerCycle=4; end
if nargin<2; Base=0; end
if nargin<1; CS=[]; IS=[]; IndVec1=[]; IndVec2=[]; return; end

[XSizeFunc,YSizeFunc,NSlice,NShot]=size(FuncMat);
NLag=MaxLag-MinLag+1;
CS=zeros(XSizeFunc,YSizeFunc,NSlice);
IS=CS;
if NLag>0
	CSMat=zeros(XSizeFunc,YSizeFunc,NSlice,NLag);
	ISMat=CSMat;
end
Message=char(zeros(6,100));
hMsgbox=msgbox(Message,'Progress Info:','help'); 
hChildMsgbox=get(hMsgbox,'Children');
set(hChildMsgbox(length(hChildMsgbox)),'Visible','off');	

ActLag=0;
for LagIndex=MinLag:MaxLag
	ActLag=ActLag+1;
	FigTitle='Mutual Information Shot Order';
	[IndVec1,IndVec2,Num,LastShot]=CalcInd12Vec(NShot,NShotPerCycle,Base+LagIndex,FigTitle);
	FuncMat1=squeeze(FuncMat(:,:,:,IndVec1));
	FuncMat2=squeeze(FuncMat(:,:,:,IndVec2));
	B=ones(Num,1);
	for SliceInd=1:NSlice
		Message=char('');
		Message=char(Message,['Calculate Mutual Information Lag ',int2str(LagIndex)]);
		Message=char(Message,'');
		Message=char(Message,['Lag   ',int2str(ActLag),' of ',int2str(NLag)]);
		Message=char(Message,['Slice ' int2str(SliceInd) ' of ' int2str(NSlice) ' ...']);
		set(hChildMsgbox(2),'String',Message);
		pause(.001);
		for i=1:XSizeFunc
			for j=1:YSizeFunc
				Tmp1=squeeze(FuncMat1(i,j,SliceInd,:))';
				Tmp2=squeeze(FuncMat2(i,j,SliceInd,:))';
				Bx1=B*Tmp1;
				Bx2=B*Tmp2;
				Dist1=(Bx1-Bx1').^2;
				Sig1=sum(sum(sqrt(Dist1)))./(Num.*(Num-1));
				Dist2=(Bx2-Bx2').^2;
				Sig2=sum(sum(sqrt(Dist2)))./(Num.*(Num-1));
				if Sig1^2>0 & Sig2^2>0
					Ker1=exp(-Dist1/(2.*Sig1^2));
					Ker2=exp(-Dist2/(2.*Sig2^2));
					DistC=(Bx2-Bx1').^2;
					KerC=exp(-DistC/(Sig1^2+Sig2^2));
					Num2=Num.^2;
					Ip1=1./Num2*sum(sum(Ker1));
					Ip2=1./Num2*sum(sum(Ker2));
					Cip=1./Num2*sum(sum(KerC));
					CS(i,j,SliceInd)=1-Cip*Cip/(Ip1*Ip2);
					IS(i,j,SliceInd)=Ip1+Ip2-2.*Cip;
				end
			end
		end
	end
	if NLag>1
		CSMat(:,:,:,ActLag)=CS;
		ISMat(:,:,:,ActLag)=IS;
	end
end
close(hMsgbox);
if NLag>1
	[CS1,CS2]=max(CSMat,[],4);
	CS2=CS2-1;
	[IS1,IS2]=max(ISMat,[],4);
	IS2=IS2-1;
else
	CS1=CS;
	CS2=[];
end
return;





