Non-negative matrix factorization

// the rows of matrix H contain the estimated components
// the output matrix NMFcomponents maps the proportions of components
function NMF(wave w3D, variable components)
    variable rows = DimSize(w3D, 0)
    variable cols = DimSize(w3D, 1)
    variable layers = DimSize(w3D, 2)
    variable pixels = rows * cols
    variable mincycles = 1000, maxcycles = 2000
    
    Duplicate/free w3D w_input
    Redimension/E=1/N=(pixels, layers) w_input
    
    w_input = max(0, w_input)
    
    // rows of W are pixels of w3D, columns are z-dimension of w3D
    Make/O/N=(pixels, components)/Y=(WaveType(w_input)) W
    Make/O/N=(components, layers)/Y=(WaveType(w_input)) H
    W = abs(gnoise(1))
    H = abs(gnoise(1))
    
    int i
    variable err_old = Inf
    variable start = ticks
    
    for (i=0;i<maxcycles;i++)
                
        MatrixOP/O W = W * ((w_input x H^t) / ((W x H) x H^t))
        MatrixOP/O H = H * ((W^t x w_input) / ((W^t x W) x H))
        
        if ((ticks-start)/60>10)
            Make/N=(rows, cols, components)/O NMFcomponents
            NMFcomponents = W[p + q * rows][r]
            DoUpdate
            DoAlert 1, num2str(i) + " iterations completed.\r\rContinue?"
            if (v_flag == 2)
                break
            endif
            start = ticks
        endif
                
        if (i >= mincycles)
            MatrixOP/O err = sqrt(sum(sq(w_input-(W x H))))
            if(err[0][0]/err_old > 0.9995)
                break
            endif
            err_old = err[0][0]
        endif
    endfor
    
    Make/N=(rows, cols, components)/O NMFcomponents
    NMFcomponents = W[p + q * rows][r]
    
    return i
end