I am new to both Matlab and Python and am converting some Matlab codes to it's Python equivalent. The issue I am facing is with converting from sparse(i,j,v,m,n) to csr_matrix((data, (row_ind, col_ind)), [shape=(M, N)]).
In this code, i, j and row_in, col_ind will be passed with an index array - idx of size(124416, 1), while v and data will be passed with a 2D array - D22 of size(290, 434)
Matlab:
...
H = 288;
W = 432;
N = (H+2)*(W+2);
mask = zeros(H+2, W+2);
mask(2:end-1, 2:end-1) = 1;
idx = find(mask==1);
>>>idx = [292, ..., 579, 582 ..., 869, ... , 125282, ..., 125569]
A = sparse(idx, idx+1, -D22(idx), N, N);
B = sparse(idx, idx-1, -D22(idx), N, N);
C = sparse(idx, idx+H+2, -D22(idx-1), N, N);
D = sparse(idx, idx-H-2, -D22(idx-1), N, N);
...
spy(A) first entry is m(293, 292) - (idx,idx+1), which was what I expected.
spy(B) m(292, 293) - (idx,idx-1). I was expecting it to be m(291, 292), believing that idx-1 would return an array [291, ..., 578, 581 ..., 868, ... , 125281, ..., 125568]
spy(C) - m(582, 292) - (idx,idx+H+2)
spy(D) - m(292, 582) - (idx,idx-H-2)
Hence, given that was how I understood the indexing order, I translated the code into Python in this form
Python:
...
H = 288
W = 432
N = (H+2) * (W+2)
mask = np.zeros([H+2, W+2])
mask[1:-1,1:-1] = 1
idx = np.nonzero(mask.transpose() == 1)
idx = np.vstack((idx[1], idx[0]))
idx = np.ravel_multi_index(idx, ((H+2),(W+2)), order='F').copy() # Linear Indexing as per Matlab
>>> idx
array([291, ..., 578, 581 ..., 868, ... , 125281, ..., 125568])
idx_ = np.unravel_index(idx, ((H+2),(W+2)), order='F') # *** Back to Linear Indexing
idx_ = np.column_stack((idx_[0], idx_[1])) # *** combine tuple of 2 arrays
idx_H_2 = np.unravel_index(idx-H-2, ((H+2),(W+2)), order='F')
idx_H_2 = np.column_stack((idx_H_2[0], idx_H_2[1]))
A = sp.csr_matrix((-D22[idx_[:,0], idx_[:,1]], (idx+1,idx)), shape = (N,N))
B = sp.csr_matrix((-D22[idx_[:,0], idx_[:,1]], (idx-1,idx)), shape = (N,N))
C = sp.csr_matrix((-D11[idx_[:,0], idx_[:,1]], (idx+H+2,idx)), shape = (N,N))
D = sp.csr_matrix((-D11[idx_H_2[:,0], idx_H_2[:,1]], (idx-H-2,idx)), shape = (N,N))
...
For A, the first entry is p(292, 291) - (idx+1,idx), and since Python starts from zero index, it refers to Matlab m(293, 292).
However for B, the first entry is p(290, 291) - (idx-1,idx), which was what I had expected (the equivalent in Matlab should be m(291, 292) ), but as mentioned earlier the Matlab code returns (292, 293) instead.
C - p(581, 291) - (idx+H+2,idx)
D - p(1, 291) - (idx-H-2,idx)
Could anyone kindly explain what I might have understood incorrectly, and how should I revise my Python code to reflect the Matlab code more accurately.
Oh and just one more qns :)
Matlab:
A = A(idx,idx);
Python:
A = A[idx,:][:,idx]
Are the equivalent?
Thank you so much for your all kind help and time.
These lines are confusing:
py(A) first entry is m(293, 292) - (idx,idx+1), which was what I expected.
spy(B) m(292, 293) - (idx,idx-1). I was expecting it to be m(291, 292), believing that idx-1 would return an array [291, ..., 578, 581 ..., 868, ... , 125281, ..., 125568]
spy(C) - m(582, 292) - (idx,idx+H+2)
spy(D) - m(292, 582) - (idx,idx-H-2)
What is m(293,292)? Why the reverse in coordinates? Is that because of how spy plots the axes? p(...) for numpy code is equally confusing. In my (smaller) samples, A, B etc all have nonzeros where I expect.
By the way, are there any zeros in D22(idx)?
In any case, you've created 4 sparse matrices, with values along one diagonal or other, with periodic zero gaps.
A(idx, idx+1) has the same nonzero values as A, but contiguously on the main diagonal.
Condensed version of the numpy code is:
In [159]: idx=np.where(mask.ravel()==1)[0]
In [160]: A=sparse.csr_matrix((np.ones_like(idx),(idx,idx+1)),shape=(N,N))
I'm ignoring the F v C order, and the D22 array. If I had D22 matrix, I'd try to use D22.ravel[idx] (to match how I create and index mask). I don't think those details matter when comparing the overall generation of the matrices and their indexing.
A.tocoo().row and A.tocoo().col is a handy way of seeing the row and column indexes of the nonzero elements. A.nonzero() does this as well (with virtually the same code).
Yes, A[idx,:][:,idx+1] produces the same submatrix.
A[idx, idx+1] gives a 1d vector of those diagonal values.
You need to transform the first index array into a 'column' vector to select a block (as the MATLAB version does):
A[np.ix_(idx,idx+1)] # or with
A[idx[:,None],idx+1]