Take the 2-minute tour ×
Code Review Stack Exchange is a question and answer site for peer programmer code reviews. It's 100% free, no registration required.

I am very much interested in the Reverse Cuthil McKee Algorithm. I have seen Fortran and C or C++ implementations of it, and I decided that it would be a nice exercise to implement it in Python. I know this algorithm is quite domain specific, but I would still be happy to see what kind of comments I get regarding:

  • Correctness - I am not sure my only test case works for others, although I did some comparison to the Octave and Matlab version.
  • Speed - Of course a C version would be faster. However, is there some Python improvements which can be done?
  • Readability - Is this code clear enough to other peer programmers?

The code:

import numpy as np

def getDegree(Graph):
    """
    find the degree of each node. That is the number
    of neighbours or connections.
    (number of non-zero elements) in each row minus 1.
    Graph is a Cubic Matrix.
    """
    degree = [0]*Graph.shape[0]
    for row in range(Graph.shape[0]):
        degree[row] = len(np.flatnonzero(Graph[row]))-1
    return degree

def getAdjcncy(Mat):
    """
    return the adjacncy matrix for each node
    """
    adj = [0]*Mat.shape[0]
    for i in xrange(Mat.shape[0]):
        q=np.flatnonzero(Mat[i])
        q=list(q)
        q.pop(q.index(i))
        adj[i] = q
    return adj

def RCM_loop(deg,start, adj,pivots,R):
    """
    Reverse Cuthil McKee ordering of an adjacency Matrix
    """
    digar=np.array(deg)
    # use np.where here to get indecies of minimums
    if start not in R:
        R.append(start)
    Q=adj[start]
    for idx, item in enumerate(Q):
        if item not in R:
            R.append(item)
    Q=adj[R[-1]]
    if set(Q).issubset(set(R)) and len(R) < len(deg) :
         p = pivots[0]
         pivots.pop(0)
         return RCM_loop(deg,p,adj,pivots,R)
    elif len(R) < len(deg):
         return RCM_loop(deg,R[-1],adj,pivots,R)
    else:
         R.reverse() 
         return R      

def test():
    """
    test the RCM loop
    """
    A = np.diag(np.ones(8))
    print A
    nzc=[[4],[2,5,7],[1,4],[6],[0,2],[1,7],[3],[1,5]]

    for i in range(len(nzc)):
        for j in nzc[i]:
            A[i,j]=1
    # define the Result queue
    R = ["C"]*A.shape[0]
    adj = getAdjcncy(A)
    degree = getDegree(A)
    digar=np.array(degree)
    pivots = list(np.where(digar == digar.min())[0])
    inl=[]
    Res = RCM_loop(degree,0, adj,pivots,inl)
    print degree
    print adj
    print "solution:", Res
    print "correct:", [6,3,7,5,1,2,4,0] 

if __name__ == '__main__':
    test()
share|improve this question
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.