1

Hello I have a 3d numpy array of shape (nc, nx, ny). I would like to plot it (with matpoltlib) using x and y as the axis, but in such a way that all c values are close to each other and arranged in a rectangle around a central value (you can assume nc is a power of 2).

Thanks to @mad-physicist's answer I think I managed to come up with a general solution, but it is not very efficient.

My code looks like this:

nc1 , nc2 = 4, 4
(nx, ny, nc) = (2,3,nc1*nc2)
data = np.reshape(np.arange(nx* ny* nc),(nc, nx, ny))

res = np.zeros((nx*nc1,ny*nc2))

for i in range(nx) :
    for j in range(ny) :
        tmp = data[:,i,j].reshape(nc1,nc2)
        res[i*nc1:(i+1)*nc1,j*nc2:(j+1)*nc2] = tmp

I am trying to find a way to avoid the loops using numpy functions. Any suggestion on how to do this?

To clarify, here is a picture of what I have in mind.

matrix transformation

muzzle
  • 315
  • 3
  • 9
  • 1
    `4 numbers correponding to x=0, y=0 (0,4,8,10)` Do you mean `(0,4,8,12)`? – Quang Hoang Jan 20 '21 at 17:54
  • Ooooops, thanks for catching this!!! I'll correct the example! – muzzle Jan 20 '21 at 17:56
  • Also, assumingly `(0,4,8,12)` has length equal to `nc`. How do you relate that with the `2x2` shape of `0,4,8,12` in the expected output in general? – Quang Hoang Jan 20 '21 at 17:57
  • I want to turn the vector of length nc into a rectangle. Assume that nc is a power of 2, so you can always turn it into a rectangle (2**i, 2**j). I think I can deal with the general case later. – muzzle Jan 20 '21 at 18:03
  • How would you do that in a loop? That's generally the best place to start, and avoids any ambiguity about what you mean. – Mad Physicist Jan 20 '21 at 18:13
  • Also, I suggest using numbers that are relatively prime for your dimensions, so it's clear where things go: `8, 3, 5` would be a good choice here. – Mad Physicist Jan 20 '21 at 18:22
  • Not sure if this is a dupe, but it may help: https://stackoverflow.com/q/62459707/2988730 – Mad Physicist Jan 20 '21 at 18:24
  • One more resquest. Some of the numbers in the original prose no longer match the code you show. Could you please update the prose? or example? whichever you think conveys the issue better. – Mad Physicist Jan 20 '21 at 22:14
  • So just to clarify, you want a 2D array that is basically a tilling of nx by ny blocks, each block of size nc1 by nc2? – Mad Physicist Jan 20 '21 at 22:20
  • I'm having trouble understanding why your output shape is `(4, 15)`, not `(6, 10)`. Please explain. I don't understand what you mean by "close" given how you rearranged the data. If you can, write out a for loop or some other recipe that explains how you got the arrangement in the final array, not just for the first pixel. – Mad Physicist Jan 20 '21 at 23:03
  • I have (hopefully) fixed my example. Sorry, @mad-physicist, the previous one had been screwed up by a cut and paste disaster. – muzzle Jan 21 '21 at 11:13
  • It's been an arduous journey, but you now have the answer you were looking for. Thanks for keeping up with all the requests for clarification. – Mad Physicist Jan 21 '21 at 15:30

2 Answers2

2

Here is a tentative approach to rearrange the memory layout closer to what you want, so you can de-interlace in a single step:

y = np.moveaxis(data.reshape(nc // 2, 2, nx, ny), 1, -1).reshape(nc, -1)
result = np.concatenate((y[::2], y[1::2]))

You may have to play with the dimensions to make sure you are rearranging based on the ones you want.

Mad Physicist
  • 107,652
  • 25
  • 181
  • 264
  • Thank you!! I think this is going in the right direction, but I cant generalize the second line to the case of nc = k * l. Do you have any ideas? – muzzle Jan 20 '21 at 20:20
  • concatenate will only work if k=2, right? In the general case the solution is more complex, but I have a hard time visualizing it. – muzzle Jan 20 '21 at 20:31
  • 1
    @muzzle. That's why I want you to make a useful MCVE – Mad Physicist Jan 20 '21 at 21:04
  • I have added a (hopefully) better example. Does that answer your questions on what transformation I am trying to achieve? – muzzle Jan 21 '21 at 11:09
1

Numpy assumes C order. So when you reshape, if your array is not C contiguous, it generally makes a copy. It also means that you can rearrange the memory layout of an array by say transposing and then reshaping it.

You want to take an array of shape (nc1 * nc2, nx, ny) into one of shape (nx * nc1, ny * nc2). One approach is to split the array into 4D first. Adding an extra dimension will not change the memory layout:

data = data.reshape(nc1, nc2, nx, ny)

Now you can move the dimensions around to get the memory layout you need. This will change your strides without copying memory, since the array is still C contiguous in the last step:

data = data.transpose(2, 0, 3, 1) # nx, nc1, ny, nc2

Notice that we moved nc1 after nx and nc2 after ny. In C order, we want the groups to be contiguous. Later dimensions are close together than earlier ones. So raveling [[1, 2], [3, 4]] produces [1, 2, 3, 4], while raveling [[1, 3], [2, 4]] produces [1, 3, 2, 4].

Now you can get the final array. Since the data is no longer C contiguous after rearranging the dimensions, the following reshape will make the copy that actually rearranges the data:

data = data.reshape(nx * nc1, ny * nc2)

You can write this in one line:

result = data.reshape(nc1, nc2, nx, ny).transpose(2, 0, 3, 1).reshape(nx * nc1, ny * nc2)

In principle, this is as efficient as you can get: the first reshape changes the shape and strides; the transpose swaps around the order of the shape and strides. Only a single copy operation happens during the last reshape.

Mad Physicist
  • 107,652
  • 25
  • 181
  • 264