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