Generating GRF with Odd vs. Even Number of Pixels

[bmbuncher]

Hi! I tried to generate 2D GRFs (using powerbox.Powerbox as shown here), but have run into issues depending on whether I use an odd or even number of pixels when generating the field. When using an even number of pixels, delta_x comes out fine; however, when using an odd number of pixels, delta_x has large top-left to bottom-right stripes. I presume that it has to do with how pyfftw handles odd vs. even numbers of modes, but it's pretty odd.

I attached a MWE, along with images of delta_x and the sampeld particles. These used P(k) ~ k^-1.2, a side length of 1, 512 (or 513) pixels, and ~10,000 particles each. I used two separate amplitudes: a "reasonable" one (where there is little clipping of the field) and an exaggerated one (which requies ~30% of the pixels in the field to be clipped to -1.0); the exaggerated one was used to show that it affects the discrete sampling.

Any suggestions (other than just using an even number of pixels)?

I'm using powerbox version v0.6.1, numpy v1.15.4, and pyfftw v0.12.0 on an Ubuntu 20.04 OS.

[bmbuncher]

Whoops! I forgot the files. Here they are. pbox_odd_even.zip

[steven-murray]

Thanks for reporting this @bmbuncher. This certainly seems like something that could happen. I'll look into it, as I don't have any ideas off the top of my head.

Can you confirm if the same happens if you uninstall pyfftw?

I guess one solution may be just to emit a warning when using an odd number of pixels (at least with pyfftw). But I'll have a crack at solving it first! If in the meantime you solve it, PRs are most welcome :-)

[bmbuncher]

I can confirm that it still has the same issue without pyfft installed

[steven-murray]

Thanks @bmbuncher, that is helpful

[bmbuncher]

I did some digging, and it appears that the issue stems from how the k-matrices are made Hermitian. If you change this so that self._even is always True, the issue vanishes (at least visually; I haven't checked how it affects things more carefully); similarly, if you set it to always be False, the anisotropy issue appears for both odd and even modes. If I understand correctly, when you have an even number of pixels, you generate a Hermitian random field array of size Npix + 1 x Npix + 1, then trim it so that you end up with an Npix x Npix array again; however, this is not done for an odd number of pixels. As a result, for an even number of pixels, you trim the row and column of your random field corresponding with the k = 0 modes; however, the fields with an odd number of pixels do retain the component that corresponds with the k = 0 mode. I ran into something similar in the past, and found that keeping the zero-mode led to issues related to the symmetries of the k-vectors, though can't give much more of an explanation until I find that piece of code.

With that said, I don't know if the solution is as simple as always removing the mode of the GRF that corresponds with k = 0, as it may affect things downstream. Secondly, there may be some other issue: with an odd number of pixels, delta_x does not have a mean of zero!

[steven-murray]

Hmmm, it's been a while since I went through the math of this code, but it seems to me that the mode that is being trimmed in the "even" case is the last mode (which in the convention of powerbox is the most positive mode). I'm pretty sure that's the intention anyway! I can't think of why that would be causing an issue, but I may not be digging deep enough. I may not have time to look at this more closely til Friday.

[bmbuncher]

That may be the case; I thought I had seen somewhere in which the field array was flipped, though I very well could be misremembering/misinterpreting it. Regardless, if it helps find the issue, I can confirm that, after performing the trimming in the line I referenced earlier, you only get a Hermitian output if you then also remove the first column and row (which, thinking about it now, would make a lot of sense if the trimming performed in the code is removing the largest positive mode), and that doing the same trimming for the negative modes yields the same result. Similarly, not doing the trimming for either also causes the same weird effect for both the odd- and even-pixel fields. This occurs whether or not pyfftw is installed.

That sounds good to me. I didn't bother looking through the math, so I'm sure that there's a reason for it. I suspect that you may know a little more about the code and math behind it than I do. If there is anything I can do in the meantime, just let me know!

[steven-murray]

Hi @bmbuncher, sorry it's taken so long to get to this, but I've done a bit more thinking, and wanted to run things by you because I'm a little confused. I think we can think of this in 1D and all the relevant concepts remain true. So what happens is that for an EVEN number of cells, we bump it up to an ODD number then produce a hermitian array, eg.

[a + bi, c + di, e + 0i, c - di, a - bi]

where the corresponding frequencies are assumed to be in purely ascending order:

[-n/2, -n/2 + 1, 0, n/2-1, n/2]

The IFFT of such an array is always purely real, as we desire. However, since the output needs to have one less cell (we put in N=4 but produced the hermitian array with n=5), we chop off the last one (which corresponds essentially to what a normal call to fftfreq will give), and then do the IFFT. However, this is NOT gauranteed to give purely real results.

So, we could do the IFFT before chopping off the last value. I'm hesitant as to whether the resulting array then has the periodic quality that would be expected, but it at least should be definitely real.

The other thing that might be slightly incorrect here is that, looking at the rfft and its inverse, they suggest that the most positive frequency should be purely real because it should be hermitian to itself -- i.e. you don't get any more information with that last cell when its odd. Currently the _make_hermitian function does not enforce this, and perhaps it should.

I wanted to get your thoughts before I dive in and change these things!

[bmbuncher]

Hi @steven-murray, no worries! Are you saying that, for the even case, you would change it so that the IFFT is done before chopping off the last value of the delta_k matrix? If so, I'm a bit confused: the code works fine with an even number of pixels, but fails with an odd number of pixels. It seems to me that chopping off the last value before taking the IFFT in the even case is what is making the code work in the first place (though it's been a while since I looked at the code thoroughly).