Simple ICA

#pragma rtGlobals=3			// Use modern global access method.
#pragma IgorVersion=6.2	// Requires Igor 6.2
	
// 
//  simpleICA(inX,reqComponents,w_init)
//  Parameters:
// 	inX is a 2D wave where columns contain a finite mix of independent components.
// 	reqComponents is the number of independent components that you want to extract.
//		This number must be less than or equal to the number of columns of inX.
//	w_init is a 2D wave of dimensions (reqComponents x reqComponents) and contains 
//  		an estimate of the mixing matrix.  You can simply pass $"" for this parameter
//		so the algorithm will use an equivalent size matrix filled with enoise().
//
//	The results of the function are the waves ICARes and WunMixing.  ICARes is a 2D 
//	wave in which each column contains an independent component.  WunMixing is a 2D
// 	wave that can be used to multiply the (re)conditioned input in order to obtain the unmixed 
//	components.
//
//	The code below implements the "deflation" approach for fastICA.  It is based on the 
//	fastICA algorithm: Hyvärinen,A (1999). Fast and Robust Fixed-Point Algorithms 
//	for Independent Component Analysis. IEEE Transactions on Neural Networks, 10(3),626-634.
// 	
//	 See testing example below the main function.
Function simpleICA(inX,reqComponents,w_init)
	Wave inX,w_init
	Variable reqComponents
	
	// The following 3 variables can be converted into function arguments.
	Variable maxNumIterations=1000
	Variable tolerance=1e-5
	Variable alpha=1
	
	Variable i,ii
	Variable iteration 
	Variable nRows=DimSize(inX,0)
	Variable nCols=DimSize(inX,1)
	
	// check the number of requested components:
	if(reqComponents>min(dimSize(inX,0),dimSize(inX,1)))
		doAlert 0,"Bad requested number of components"
		return 0
	endif
	
	// Never mess up the original data
	Duplicate/O/Free inX,xx					

	// Initialize the w matrix if it is not provided.	
	if(WaveExists(w_init)==0)
		Make/O/N=(reqComponents,reqComponents) w_init=enoise(1)
	endif
	
	// condition and transpose the input:
	MatrixOP/O xx=(NormalizeCols(subtractMean(xx,1)))^t

	// Just like PCA:
	MatrixOP/O/Free V=(xx x (xx^t))/nRows
	MatrixSVD V
	// M_VT is not used here.
	Wave M_U,W_W,M_VT									
	W_W=1.0/sqrt(w_w)
	MatrixOP/O/Free D=diagonal(W_W)
	MatrixOP/O/FREE K=D x (M_U^t)			
	KillWaves/z W_W,M_U,M_VT			 

	Duplicate/Free/R=[0,reqComponents-1][] k,kk
	Duplicate/O/FREE kk,k									

	// X1 could be output as PCA result.	
	MatrixOP/O/FREE X1=K x xx								
	// create and initialize working W; this is not an output!
	Make/O/Free/N=(reqComponents,reqComponents) W=0						
					
	for(i=1;i<=reqComponents;i+=1)										
		MatrixOP/O/FREE lcw=row(w_init,i-1)
               // decorrelating 							
		if(i>1)													
			Duplicate/O/Free lcw,tt									
			tt=0												
			for(ii=0;ii<i;ii+=1)
				MatrixOP/O/Free r_ii=row(W,ii)				// row ii of matrix W		
				MatrixOP/O/FREE ru=sum(lcw*r_ii)			// dot product		
				Variable ks=ru[0]
				MatrixOP/O/Free tt=tt+ks*r_ii							
			endfor
			MatrixOP/O/FREE lcw=lcw-tt								
		endif
		MatrixOP/O/Free lcw=normalize(lcw)						
		// iterate till convergence:	
		for(iteration=1;iteration<maxNumIterations;iteration+=1)				
			MatrixOP/O/Free wxProduct=lcw x x1						
			// should be supported by matrixop :(
			Duplicate/O/Free wxProduct,gwx
			gwx=tanh(alpha*wxProduct)									
			Duplicate/Free/R=[reqComponents,nRows] gwx,gwxf				 
			Make/O/Free/N=(reqComponents,nRows) gwxf
			// repeat the values from the first row on.
			gwxf=gwx[q]										
			Duplicate/O/FREE gwxf,gwx
			MatrixOP/O/Free x1gwxProd=x1*gwx							 
			Duplicate/O/FREE wxProduct,gwx2									 
			gwx2=alpha*(1-(tanh(alpha*wxProduct))^2)
			Variable theMean=mean(gwx2)
			MatrixOP/O/Free    wPlus=(sumRows(x1gwxProd)/numCols(x1gwxProd))^t-theMean*lcw	
			// reduce components						 
			Redimension/N=(1,reqComponents) wPlus				
			// starting from the second component;
			if(i>1)												
				Duplicate/O/FREE wPlus,tt									 
				tt=0												 
				for(ii=0;ii<i;ii+=1)					                  
					MatrixOP/O/Free r_ii=row(W,ii)				 
					MatrixOP/O/FREE ru=wPlus.(r_ii^t)							 
					ks=ru[0]
					MatrixOP/O tt=tt+ks*r_ii							 
				endfor										            
				wPlus=wPlus-tt							       			 
			endif
			MatrixOP/O/FREE wPlus=normalize(wPlus)							 
			MatrixOP/O/Free limV=abs(mag(sum(wPlus*lcw))-1)		
			printf "Iteration %d, diff=%g\r",iteration,limV[0]
			lcw=wPlus
			if(limV[0]<tolerance)
				break
			endif
		endfor
		// store the computed row in final W.
        	W[i-1][]=lcw[q]													
	endfor			// loop over components

	//  Calculate the un-mixing matrix
	MatrixOP/O WunMixing=W x K					
	// 	Un-mix; 					
	MatrixOP/O ICARes=(WunMixing x xx)^t								
End


Function testICA()
	// Create distinct frequencies:
	Make/O/N=(1000,3) ddd
	ddd[][0]=sin(2*pi*x/13)
	ddd[][1]=sin(2*pi*x/17)
	ddd[][2]=sin(2*pi*x/23)
	// Create mixing matrix:
	Make/O/N=(3,3) AA
	AA[0][0]= {0.291,0.6557,-0.5439}
	AA[0][1]= {0.5572,0.3,-0.2}
	AA[0][2]= {-0.1,-0.7,0.4}
	// Do your mixing:
	MatrixOP/O xx=ddd x AA
	// Try the ICA
	simpleICA(xx,3,$"")
	// Plot the results:
	Wave ICARes
	Display ICARes[][0]
	Display ICARes[][1]
	Display ICARes[][2]
End