There are several ways to get submatrix in numpy:

In [35]: ri = [0,2]
    ...: ci = [2,3]
    ...: a[np.reshape(ri, (-1, 1)), ci]
Out[35]: 
array([[ 2,  3],
       [10, 11]])

In [36]: a[np.ix_(ri, ci)]
Out[36]: 
array([[ 2,  3],
       [10, 11]])

In [37]: s=a[np.ix_(ri, ci)]

In [38]: np.may_share_memory(a, s)
Out[38]: False

note that the submatrix you get is a new copy, not a view of the original mat.

It probably is more useful to work with the transpose of W rather than W itself, both for human-readability and to facilitate writing the code. This means that the entries that affect each S_i are grouped together in one of the inner parentheses of W, i.e. in a row of W rather than a column as you have it now.

Then, S_i = np.array[S[j,:] for j in np.shape(S)[0] if W_T[i,j] == 1], where W_T is the transpose of W. If you need/want to stick with W as is, you need to reverse the indices i and j.

As for the outer loop, you could try to nest this in another similar comprehension without an if statement--however this might be awkward since you aren't actually building one output matrix (the S_i can easily be different dimensions, unless you're somehow guaranteed to have the same number of 1s in every column of W). This in fact raises the question of what you want--a list of these arrays S_i? Otherwise if they are separate variables as you have it written, there's no good way to refer to them in a generalizable way as they don't have indices.

Give np.ix_ a try:

Y[np.ix_([0,3],[0,3])]

This returns your desired result:

In [25]: Y = np.arange(16).reshape(4,4)
In [26]: Y[np.ix_([0,3],[0,3])]
Out[26]:
array([[ 0,  3],
       [12, 15]])

to answer your question s is only a 1d array ... (even if you did actually transpose it ... which you did not)

>>> u,s,v = linalg.svd(A)
>>> s
array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.])
>>>

for selecting a submatrix I think this does what you want ... there may be a better way

>>> rows = range(10,15)
>>> cols = range(5,8)
>>> A[rows][:,cols]
array([[ 0.,  0.,  0.],
       [ 0.,  0.,  0.],
       [ 0.,  0.,  0.],
       [ 0.,  0.,  0.],
       [ 0.,  0.,  0.]])

or probably better

>>> A[15:32, 2:7]
array([[ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.]])

With A as the input array containing those sub-matrices, you could use some reshaping and permute dimensions with np.transpose, like so -

A.reshape(2,2,3,3).transpose(0,2,1,3).reshape(-1,6)

Explanation -

• Cut the first axis into two parts : A.reshape(2,2,3,3).
• Swap the second and third axes : .transpose(0,2,1,3).
• Finally, reshape to have 6 columns : .reshape(-1,6).

Thus, the solution could be generalized like so -

m,n,r = A.shape
out = A.reshape(-1,2,n,r).transpose(0,2,1,3).reshape(-1,2*r)

Sample run

Let's re-create your A -

In [54]: A = np.arange(36).reshape(-1,3,2,3). transpose(0,2,1,3).reshape(4,3,3)

In [55]: A
Out[55]: 
array([[[ 0,  1,  2],
        [ 6,  7,  8],
        [12, 13, 14]],

       [[ 3,  4,  5],
        [ 9, 10, 11],
        [15, 16, 17]],

       [[18, 19, 20],
        [24, 25, 26],
        [30, 31, 32]],

       [[21, 22, 23],
        [27, 28, 29],
        [33, 34, 35]]])

Then, run our code -

In [56]: A.reshape(2,2,3,3).transpose(0,2,1,3).reshape(-1,6)
Out[56]: 
array([[ 0,  1,  2,  3,  4,  5],
       [ 6,  7,  8,  9, 10, 11],
       [12, 13, 14, 15, 16, 17],
       [18, 19, 20, 21, 22, 23],
       [24, 25, 26, 27, 28, 29],
       [30, 31, 32, 33, 34, 35]])

This function fills in a circle of 1s that looks pretty good.

def fill_cell(cell, corner, rad):
    m, n = cell.shape
    ctr = corner[0]+m/2, corner[1]+n/2
    x = np.arange(m) - ctr[0]
    y = np.arange(n) - ctr[1]
    X,Y = np.meshgrid(x, y, order='ij')  # could try order='xy'
    Z = ((X**2 + Y**2)<= rad**2).astype(cell.dtype)
    return Z
empty_lattice[:] = fill_cell(empty_lattice, (x,y),side/2)

Position in empty_lattice is not right - because of a difference in how you are defining the x,y coordinates and my assumptions, but I think you can sort that out.

Radius looks good, though it might be off by an integer.

To fill in multiple circles, you could iterate over the x,y values, and assign lattice values for a slice (view)

xyslice = slice(x,15), slice(y,15)
empty_lattice[xyslice] = fill_cell(empty_lattice[xyslice],...)

For overlapping circles I'd look into some sort of logical assignment

empty_lattice[xyslice] |= fill_cell(...)