Odd appearance of some photon shooting outputs in the MultiObjectDemo output

[barnabytprowe]

As discovered and discussed a little on the Pull Request #206 (for Issue #175), Mike found one bug causing difference between the DFT and photon output but (quoting Mike):

...there still seems to be a problem with the image. In fact all of the SBInterpolatedImage photon shooting images (61-100) look a bit wrong. So there might still be some bug in that code somewhere.

This Issue is to investigate this behaviour. Distinguishing things that the PSFs used in 61-100 share is that they both use SBInterpolated images, and they both use Convolve to combine optical and atmospheric parts.

[barnabytprowe]

Unfortunately I won't be able to look at this myself today, but is there agreement that it might be good to wait to see if this can be quickly resolved before generating a milestone3 tag for the code? I think so, but @rmandelb and @rmjarvis (and others) do you agree?

[rmandelb]

Yes, I agree.

[rmjarvis]

I can now more specifically quantify what my eye was telling me. (Not sure if it was obvious to everyone what I meant by "a bit wrong".) I have a script in branch '#208' called Test.py (just in the main GalSim directory -- probably don't need to keep this after we diagnose this issue).

It runs Kolmogorov * Gaussian 500 times with different noise realizations, but everything else the same. Then it looks at the statistics of the pixel outputs from those runs. There is no sky_level in this run, so the noise is purely photon noise, so it should obey Poisson statistics with var = mean. Here is what I get for the central pixel:

fft mean = 438.290000, variance = 444.174248
    range = 382.000000, 502.000000
    quartiles = 424.000000, 453.000000
phot mean = 436.862023, variance = 2458.091936
     range = 294.711212, 559.447937
     quartiles = 400.818390, 470.732635

FFT is close, but clearly photon shooting has a much larger variance than it should have.

To prove that this is specific to SBInterpolatedImage, and not a general photon shooting problem, here are the results if I switch the PSF to Moffat:

fft mean = 394.374000, variance = 395.745615
    range = 343.000000, 453.000000
    quartiles = 382.000000, 408.000000
phot mean = 394.120000, variance = 393.809218
     range = 338.000000, 453.000000
     quartiles = 381.000000, 407.000000

Right on the money.

The next step is to vary some parameters and see what exacerbates the problem and what (if anything) makes it go away. @gbernstein Gary, if you have any ideas of places to look for the bug, please let me know. It's not related to the previous bug, since it happens with dx=1 too. (This script keeps the same pixel_scale=0.28 as MultiObjectDemo.py, but I tried switching it to 1 as well, and the effect was the same.)

[rmjarvis]

Update: I'm pretty sure I found the problem. The problem was with the getPositiveFlux and getNegativeFlux for SBInerpolatedImage. They only returned the + and - flux in the image itself. But since that gets convolved with the interpolant, we need to adjust those values the same way that SBConvolve does.

I'm still tracking down some of the consequences of this fix, but I should have a solution tonight or tomorrow I think.

[rmjarvis]

I think I fixed that bug, but the variances still aren't quite right. Much closer than before, but the central pixels are now getting a bit too low a variance and the outer pixels are too high. Here is some output from my Test.py script:

pixel 0,-10:
fft mean = 0.203200, variance = 0.207551
phot mean = 0.192431, variance = 0.672496
pixel 0,-9:
fft mean = 1.046800, variance = 1.084827
phot mean = 1.055481, variance = 3.075607
pixel 0,-8:
fft mean = 5.118600, variance = 5.152765
phot mean = 5.058199, variance = 12.321863
pixel 0,-7:
fft mean = 22.736200, variance = 22.643538
phot mean = 22.895952, variance = 44.466339
pixel 0,-6:
fft mean = 87.095400, variance = 90.100319
phot mean = 87.306313, variance = 128.912483
pixel 0,-5:
fft mean = 281.602800, variance = 287.866605
phot mean = 282.328576, variance = 348.925652
pixel 0,-4:
fft mean = 749.153200, variance = 758.555441
phot mean = 748.749452, variance = 829.133601
pixel 0,-3:
fft mean = 1617.639200, variance = 1653.868997
phot mean = 1616.872439, variance = 1615.948964
pixel 0,-2:
fft mean = 2812.926800, variance = 2830.586359
phot mean = 2811.013835, variance = 2495.616676
pixel 0,-1:
fft mean = 3923.949400, variance = 3860.845009
phot mean = 3923.295635, variance = 3363.758414
pixel 0,0:
fft mean = 4384.689400, variance = 4507.176763
phot mean = 4384.671641, variance = 3941.471799
pixel 0,1:
fft mean = 3924.790800, variance = 3960.102656
phot mean = 3923.170606, variance = 3463.400559
pixel 0,2:
fft mean = 2812.836200, variance = 2776.493068
phot mean = 2812.786801, variance = 2632.982588
pixel 0,3:
fft mean = 1617.704600, variance = 1631.887716
phot mean = 1617.574637, variance = 1546.822379
pixel 0,4:
fft mean = 749.676000, variance = 746.079440
phot mean = 749.046336, variance = 823.281183
pixel 0,5:
fft mean = 281.576800, variance = 284.295361
phot mean = 281.947395, variance = 364.849849
pixel 0,6:
fft mean = 87.227400, variance = 85.120713
phot mean = 87.282819, variance = 132.350165
pixel 0,7:
fft mean = 22.636800, variance = 23.066699
phot mean = 22.726006, variance = 44.429719
pixel 0,8:
fft mean = 5.177400, variance = 5.030935
phot mean = 5.171511, variance = 12.817201
pixel 0,9:
fft mean = 1.047800, variance = 1.040923
phot mean = 1.044439, variance = 2.969668
pixel 0,10:
fft mean = 0.194400, variance = 0.196648
phot mean = 0.188831, variance = 0.676127
fft flux mean = 125996.128125, variance = 125742.036192
phot flux mean = 125991.253155, variance = 122452.752397

I'm not sure what to make of that. The means are all coming out right, just the variances are off. And the total flux (the last two lines) have the right variance, so I think this means I'm picking the overall number of photons correctly. But the variance within the image isn't quite right.

[rmjarvis]

For this new problem of the variance of the pixels not coming out right, I've made a little progress. I've discovered that it seems to be particular to Lanczos, and it gets worse for larger n values.

However, this also correlates with larger eta values (ratio of the number of negative photons to the total), so it might not be Lanczos per se -- it might just be something getting worse as eta increases. But Lanczos n=1 has eta=0 and it also seems to be wrong (although it also doesn't get the mean right, so there might be something else going on with that case). And Quintic has a relatively large eta, but it doesn't show any significant error in the variance.

Here is what I get for the central pixel for different interpolants (along with their eta): (The mean and variance should be equal.)

• Delta: eta = 0, phot mean = 438.136200, variance = 441.961642
• Nearest: eta = 0, phot mean = 433.113200, variance = 438.380062
• Linear: eta = 0, phot mean = 428.787200, variance = 431.640644
• Cubic: eta = 0.133, phot mean = 438.430343, variance = 422.354481
• Quintic: eta = 0.184, phot mean = 438.439510, variance = 433.095946
• Lanczos n=1: eta = 0, phot mean = 403.417644, variance = 379.371326
• Lanczos n=3: eta = 0.229, phot mean = 437.993522, variance = 404.570501
• Lanczos n=5: eta = 0.295, phot mean = 438.679754, variance = 395.059243
• Lanczos n=7: eta = 0.327, phot mean = 438.415681, variance = 329.288721
• Lanczos n=15: eta = 0.376, phot mean = 438.494012, variance = 268.876374

So I think there is probably something wrong with the implementation of Lanczos. I didn't delve into the details of the math, but I did tweak some of the arbitrary magic numbers in the code (like how big the uStep size is for the table, how high a max u value to use in the table, etc.) and none of them made any significant difference.

Gary, could you maybe take a look at this and see if you can think what might be going on here?

[gbernstein]

Hi Mike - I'll put a reply here for the record even though I could also walk the 20 feet to discuss this in person. I am not convinced just yet that this is a bug or a problem. I presume you're choosing the number of photons to shoot based on the formula with eta that we derived week before last. Maybe it is the case that the central pixel has less pos/neg cancellation than the average of the entire image, and therefore it has less fractional variance than the outer parts of the image and less than you have specified in your input. Would that explain it? Another possible test would be to use an SBInterpolatedImage that just has a single nonzero pixel interpolated with Lanczos5, and look at the statistics of an image obtained by photon-shooting with dx=0.25 input pixels. That might help figure out whether this is a bug or just a consequence of the somewhat unpredicatable statistical properties of photon shooting with pos & neg photons.

[rmjarvis]

Your intuition was right on this. The edges do have higher eta values than the center. Hence more variance. I guess this is because the Lanczos filter is fairly broad in its extent, so the negative parts are multiplied by the larger flux nearer the center.

However, I'm not sure that this means that it's not a problem. The statistics of the resulting image are different for photon shooting and for fft. It seems like this difference would matter for sensitive tests of shear measurements, since we know that anisotropic noise can cause biases.

Does this mean that when photon shooting things with eta > 0, we need to use many more photons than N = flux, and then add the noise in afterward? Or at least use somewhat more photons and then try to map out the eta value for each pixel and add in extra noise to make each pixel have the expected variance?

[gbernstein]

On Jul 2, 2012, at 4:12 PM, Mike Jarvis wrote:

Does this mean that when photon shooting things with eta > 0, we need to use many more photons than N = flux, and then add the noise in afterward? Or at least use somewhat more photons and then try to map out the eta value for each pixel and add in extra noise to make each pixel have the expected variance?

I do think it means that you cannot issue a guaranteed exact noise level for photon shooting unless everything is non-negative (e.g. using bilinear interpolation on a non-negative SBInterpolatedImage). I figured that the two use cases for photon shooting were:

a) check the FFT-based methods for accuracy. In this case we do not need to know exactly the noise level, we will just be shooting enough photons to bring the noise down to roughly the level desired for some test. And anyway we might be happy limiting these tests to cases that are non-negative.

b) draw images more quickly at low S/N. In this case, if there are negative photons, then yes, I would think you need to shoot enough photons to make the noise well below the target, then add in noise per pixel to bring it to the desired level. But in the case where we are sky-limited, not limited by the shot noise of the source itself, is this already what we're doing?

I think that you need to draw a bunch of photons or do an FFT to figure out the eta of each pixel. Once you've done that, you already have a high-S/N version of your image so you might as well just add the noise with a deviate instead of shooting photons again. So is there any utility to trying to fix up the somewhat unpredictable noise levels of photon-shooting?

[rmjarvis]

OK. I guess these comments are relevant to our recommendations to users about when and how to use photon shooting. cf. Issue #204.

The other implication is that I think this means we should migrate AtmosphericPSF and OpticalPSF into C++ where they don't need to be implemented using SBInterpolatedImage. They are currently implemented as fourier transforms of a sampled image in k-space, which I think means that sinc interpolation (or Lanczos for a good approximation thereof) is the right interpolant. But if we implement it directly in C++, then we don't need to do this.

[barnabytprowe]

Pull request merged, closing...