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?
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 = *Graph.shape for row in range(Graph.shape): degree[row] = len(np.flatnonzero(Graph[row]))-1 return degree def getAdjcncy(Mat): """ return the adjacncy matrix for each node """ adj = *Mat.shape for i in xrange(Mat.shape): 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 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=[,[2,5,7],[1,4],,[0,2],[1,7],,[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 adj = getAdjcncy(A) degree = getDegree(A) digar=np.array(degree) pivots = list(np.where(digar == digar.min())) 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()