%==================================================================
%
%	PCAICAMusic
%
%	Andreas Keil & Markus Junghoefer	[1998]
%   This software is protected by german copyright and international treaties.             
%   Copyright 2004 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.                                                                      
%
% 
%	Function definition
%
	
	function PCAICAMusic(action);

%
%
%=================================================================
global hPlot2dMenuFig VMat UMat
hList=get(hPlot2dMenuFig,'Userdata');
Domain=get(hList(20),'Value');
ActFileNo=get(hList(1),'Userdata');
NChanVec=get(hList(16),'Userdata');
NChan=NChanVec(1);
FileNoVec=get(hList(56),'Userdata');
NListHdr=get(hList(51),'Userdata');
MinPoint=get(hList(54),'Value');
MaxPoint=get(hList(55),'Value');
NPoints=MaxPoint-MinPoint+1;
TrigPoint=get(hList(22),'Value');
SampRate = get(hList(60),'Userdata');
MinPointMs = (MinPoint-TrigPoint).*1000./SampRate;	
MaxPointMs = (MaxPoint-TrigPoint).*1000./SampRate;
CalcType=[action(length(action)-11:length(action)-9)];
AllActType=[action(length(action)-8:length(action)-6)];
OrgActType=[action(length(action)-5:length(action)-3)];
DataType=[action(length(action)-2:length(action))];
if strcmp(CalcType,'PCA')
	if strcmp(DataType,'Cov')
		VFigureName='Time Courses of Covariance Matrix PCA';
		SFigureName='Variance Vector of Covariance Matrix PCA';
	elseif strcmp(DataType,'Cor')
		VFigureName='Time Courses of Correlation Matrix PCA';
		SFigureName='Variance Vector of Correlation Matrix PCA';
	else
		VFigureName='Time Courses of Data Matrix PCA';
		SFigureName='Variance Vector of Data Matrix PCA';
	end
elseif strcmp(CalcType,'ICA')
	if strcmp(DataType,'Cov')
		VFigureName='Time Courses of Covariance Matrix ICA';
		SFigureName='Variance Vector of Covariance Matrix ICA';
	elseif strcmp(DataType,'Cor')
		VFigureName='Time Courses of Correlation Matrix ICA';
		SFigureName='Variance Vector of Correlation Matrix ICA';
	else
		VFigureName='Time Courses of Data Matrix ICA';
		SFigureName='Variance Vector of Data Matrix ICA';
	end
elseif strcmp(CalcType,'VMA')
	if strcmp(DataType,'Cov')
		VFigureName='Time Courses of VariMax of Covariance Matrix';
		SFigureName='Variance Vector of VariMax of Covariance Matrix';
	elseif strcmp(DataType,'Cor')
		VFigureName='Time Courses of VariMax of Correlation Matrix';
		SFigureName='Variance Vector of VariMax of Correlation Matrix';
	else
		VFigureName='Time Courses of VariMax of Data Matrix';
		SFigureName='Variance Vector of VariMax of Data Matrix';
	end
elseif strcmp(CalcType,'PMA')
	if strcmp(DataType,'Cov')
		VFigureName='Time Courses of ProMax of Covariance Matrix';
		SFigureName='Variance Vector of ProMax of Covariance Matrix';
	elseif strcmp(DataType,'Cor')
		VFigureName='Time Courses of ProMax of Correlation Matrix';
		SFigureName='Variance Vector of ProMax of Correlation Matrix';
	else
		VFigureName='Time Courses of ProMax of Data Matrix';
		SFigureName='Variance Vector of ProMax of Data Matrix';
	end
end
if strcmp(AllActType,'Act');
	FileNoVec=FileNoVec(ActFileNo);
end

if isempty(MinPointMs); MinPointMs=1; end
if isempty(MaxPointMs); MaxPointMs=MinPointMs+NPoints-1; end
if strcmp(computer,'MAC2')	
   VFigPosVec=[0.4 0 0.6 1];	
	SFigPosVec=[0.6 0 0.4 0.4];
elseif strcmp(computer,'PCWIN') 
  	VFigPosVec=[0.405 0.05 .60 0.88];  
  	SFigPosVec=[0.61 0.05 0.4 0.4];
elseif strcmp(computer,'SOL2')
	VFigPosVec=[0.41 0.01 .601 0.935];  
  	SFigPosVec=[0.6224 0.01 0.4 0.4];		
end;
NComp1=1;
for	FileIndex=1:length(FileNoVec)
	if ~strcmp(AllActType,'Cel')
		if strcmp(OrgActType,'Act'); 
			ActTmpMat=get(hList(NListHdr+(FileNoVec(FileIndex)-1).*50+1),'Userdata');
		elseif strcmp(OrgActType,'Org'); 
			ActTmpMat=get(hList(NListHdr+(FileNoVec(FileIndex)-1).*50+24),'Userdata');
		end
		[NChan,NPlotPoints]=size(ActTmpMat);
	else
		FilePath=get(hList(NListHdr+(ActFileNo-1).*50+4),'Userdata');
		EgisAvgVec=get(hList(NListHdr+(ActFileNo-1).*50+31),'Userdata');
		NCells=EgisAvgVec(3);
		ActTmpMat=zeros(NChan,NPoints.*NCells);
		[ActTmpMat]=ReadEGISAvgFile(FilePath,1,0);	%Read all cells in a row => [NChan,NPoints*NCells]=size(ActTmpMat)
		[NChan,NPoints]=size(ActTmpMat);
		MaxPoint=NPoints;
		NPlotPoints=NPoints;
	end
	if FileIndex==1
   		MinPointsChan=min(NPoints,NChan);
		UMat=zeros(NChan,NChan,length(FileNoVec));
		if strcmp(CalcType,'PCA')
			VMat=zeros(NPoints,NPoints,length(FileNoVec));
		else
			VMat=zeros(NChan,NPoints,length(FileNoVec));
		end
		VarVecMat=zeros(MinPointsChan,length(FileNoVec));
  	end
	ActTmpMat=ActTmpMat(:,MinPoint:MaxPoint);
	if strcmp(DataType,'Cov')
		ActTmpMat = cov(ActTmpMat);
	elseif strcmp(DataType,'Cor')
		ActTmpMat = corrcoef(ActTmpMat);	
	end
	if strcmp(CalcType,'PCA')
		[U,S,V] = svd(ActTmpMat);
		V=V';
		VarVec = (diag(S)).^2/sum(diag(S).^2).*100;
		NComp1=max(2,1+length(find(VarVec>=1)));
	elseif strcmp(CalcType,'PMA')
		NComp1=4
		ProMaxRotMat=promax(ActTmpMat',NComp1,0,10^-10);
		RotData=ProMaxRotMat*ActTmpMat';        % compute rotated data
		size(RotData)
		U=zeros(NChan,NChan);
		U(:,1:NComp1)=RotData';
		V=zeros(NChan,NPoints);
		V(1:NComp1,:)=(RotData*ActTmpMat)./NChan;
		VarVec=ones(MinPointsChan,1);
	elseif strcmp(CalcType,'VMA')
		NComp1=4
		[eigenvectors,eigenvalues,compressed,datamean] = pcsquash(ActTmpMat',NComp1);
		[VariMaxRotMat,RotData] = varimax(compressed,10^-20); % run Varimax on the (projected) data
		U=zeros(NChan,NChan);
		U(:,1:NComp1)=RotData';
		V=zeros(NChan,NPoints);
		V(1:NComp1,:)=(RotData*ActTmpMat)./NChan;
		VarVec=ones(MinPointsChan,1);
	elseif strcmp(CalcType,'ICA')
% 		ActTmpMat=ActTmpMat(1:128,:);
% 		fprintf(1,'Actmp=Actmp(1:128,:);')
		[Weights,Sphere]=runica(ActTmpMat);
		UM=Weights*Sphere; %UnmixingMatrix = Weights*Sphere {sphering off -> eye(chans)
		if any(~isreal(UM))
	 		V=abs((UM*ActTmpMat));
			U=abs(inv(UM)); 
		else
			V=(UM*ActTmpMat);
			U=inv(UM); 
		end
		VarVec=(diag(U)).^2/sum(diag(U).^2);
		NComp1=MinPointsChan;
	end
	UMat(:,:,FileIndex)=U;
	VMat(:,:,FileIndex)=V;
    VarVecMat(:,FileIndex)=VarVec;
	MinVVec(FileIndex)=min(min(V(1:NComp1,:)));
   	MaxVVec(FileIndex)=max(max(V(1:NComp1,:)));
end
set(hList(113),'Userdata',VarVecMat);
[SigValVec] = FindRCrit(MinPointsChan);
set(hList(114),'Userdata',SigValVec);
if NComp1==1
	NRow1=1; NCol1=1;
elseif NComp1==2
	NRow1=1; NCol1=2;
elseif NComp1==3 | NComp1==4
	NRow1=2; NCol1=2;
elseif NComp1==5 | NComp1==6
	NRow1=2; NCol1=3;
else
	NRow1=ceil(sqrt(NComp1));
	NCol1=NRow1;
end
hPCAVFig=figure('Units','normal','Position',VFigPosVec, ...
       	'Name',VFigureName,'Color',[1 1 1],'NumberTitle','off','Visible','on');
TimeVec=linspace(MinPointMs,MaxPointMs,NPoints);
for Comp = 1 : NComp1
	hPCAVhAxesVec(Comp)=subplot(NRow1,NCol1,Comp);
end

hPCASFig=figure('Units','normal','Position',SFigPosVec, ...
       	'Name',SFigureName,'Color',[1 1 1],'NumberTitle','off','Visible','on');
if MinPointsChan>NComp1
	hPCAShAxesVec(1)=subplot(1,3,1);
	hPCAShAxesVec(2)=subplot(1,3,2);
	hPCAShAxesVec(3)=subplot(1,3,3);
else
	hPCAShAxesVec(1)=subplot(1,2,1);
	hPCAShAxesVec(3)=subplot(1,2,2);
end
SumVarVecMat=zeros(size(VarVecMat));
for	FileIndex=1:length(FileNoVec)
  	LineStyleStr=get(hList(NListHdr+(FileNoVec(FileIndex)-1).*50+17),'Userdata');
	ColorVal=get(hList(NListHdr+(FileNoVec(FileIndex)-1).*50+16),'Value');
	if 		ColorVal==1; 	ColorStr='k'; 	
	elseif 	ColorVal==2; 	ColorStr='w'; 	
   	elseif 	ColorVal==3;	ColorStr='y'; 	
   	elseif 	ColorVal==4; 	ColorStr='r'; 	
	elseif 	ColorVal==5; 	ColorStr='g'; 	
	elseif 	ColorVal==6; 	ColorStr='b'; 	
	elseif 	ColorVal==7; 	ColorStr='m'; 	
	elseif 	ColorVal==8; 	ColorStr='c'; 	
	else;					ColorStr='b'; 	
	end
  	ColorLineStyleStr=[ColorStr,LineStyleStr];
	ColorStarStyleStr=[ColorStr '*'];
	VarVecVec=VarVecMat(:,FileIndex);
	if any(VarVecVec-ones(size(VarVecVec)))
		[SortVarVecVal,SortVarVecInd]=sort(VarVecVec);
		SortVarVecVal=SortVarVecVal(MinPointsChan-NComp1+1:MinPointsChan);
		SortVarVecInd=SortVarVecInd(MinPointsChan-NComp1+1:MinPointsChan);
	else
		SortVarVecInd=[NComp1:-1:1];
		SortVarVecVal=ones(NComp1,1);
	end
  	for Comp=NComp1:-1:1
		axes(hPCAVhAxesVec(NComp1-Comp+1))
		plot(TimeVec,VMat(SortVarVecInd(Comp),:,FileIndex),ColorLineStyleStr)
		hold on
		plot([min(TimeVec) max(TimeVec)],[0 0],'k--')
		if MinPointMs<0
			plot([0 0],[min(TimeVec) max(TimeVec)],'k--');
		end
    	TitleString =['Var.: ',num2str(round(SortVarVecVal(Comp).*10)./10),'% ;  Comp. # ',num2str(SortVarVecInd(Comp))];
    	xlabel('Time [ms]') 
		text('Units','Normal','FontSize',[10],'HorizontalAlignment','Center','Position',[.5 1+(length(FileNoVec)-FileIndex+1).*.04],'Color',ColorStr,'String',TitleString);
    	axis([MinPointMs MaxPointMs min(MinVVec) max(MaxVVec)])  
		SumVarVecMat(NComp1-Comp+1,FileIndex)=sum(SortVarVecVal(NComp1:-1:Comp));
   	end
	axes(hPCAShAxesVec(1))
	plot(VarVecVec,ColorLineStyleStr)
	hold on
	plot(VarVecVec,ColorStarStyleStr)
	axis([1 MinPointsChan 0 1.1.*max(VarVecMat(:,FileIndex))])
	xlabel('Component []') 
	ylabel('Variance [%]') 
	if MinPointsChan>NComp1
		axes(hPCAShAxesVec(2))		
		plot(VarVecVec,ColorLineStyleStr)
		hold on
		plot(VarVecVec,ColorStarStyleStr)
		axis([1 NComp1 0 1.1.*max(VarVecMat(:,FileIndex))])
		xlabel('Component #')
		ylabel('Variance [%]')
	end
	axes(hPCAShAxesVec(3));		
	plot(SumVarVecMat(:,FileIndex),ColorLineStyleStr)
	hold on
	plot(SumVarVecMat(:,FileIndex),ColorStarStyleStr)
	axis([1 NComp1 0 1.1.*max(SumVarVecMat(:,FileIndex))])
	xlabel('Component #') 
	ylabel('Sum of Variance [%]') 
	if strcmp(DataType,'Dat')
		ActTmpMat=zeros(NChan,NPlotPoints);
		for CompInd=1:MinPointsChan
			ActTmpMat(:,CompInd)=UMat(:,CompInd,FileIndex).*VarVecMat(CompInd,FileIndex);
		end
		set(hList(NListHdr+(FileNoVec(FileIndex)-1).*50+1),'Userdata',ActTmpMat);
		MinTmpVec(FileIndex)=min(min(ActTmpMat(:,1:MinPointsChan)));
		MaxTmpVec(FileIndex)=max(max(ActTmpMat(:,1:MinPointsChan)));
	end
end
hold off
if strcmp(DataType,'Dat')
	MinAmp=1.1.*min(MinTmpVec);
	MaxAmp=1.1.*max(MaxTmpVec);
	MinAmpVec=get(hList(52),'Userdata');
	MaxAmpVec=get(hList(53),'Userdata');
	MinPointVec=get(hList(54),'Userdata');
	MaxPointVec=get(hList(55),'Userdata');
	MinAmpVec(5)=MinAmp;
	MaxAmpVec(5)=MaxAmp;
	MinPointVec(5)=1;
	MaxPointVec(5)=NComp1;
	set(hList(52),'Value',MinAmp,'Userdata',MinAmpVec,'String',sprintf('%g',MinAmp)); 	%MinAmp
	set(hList(53),'Value',MaxAmp,'Userdata',MaxAmpVec,'String',sprintf('%g',MaxAmp)); 	%MaxAmp
	set(hList(54),'Value',1,'Userdata',MinPointVec,'String','1'); 						%MinPoint
	set(hList(55),'Value',NComp1,'Userdata',MaxPointVec,'String',sprintf('%g',NComp1)); %MaxPoint 
	if Domain~=5
		set(hList(20),'Value',5);	%Domain=5; PCA
		Plot2d('setPCADomain')
	end
	Plot2d('PlotData')
end

return;
