numerical issue for some realizations of atmosphere

[esheldon]

I'm seeing large offsets in the photon x, y values, as well as the occasional nan, for some realizations of an atmosphere. The offsets can be very high, like 10 million pixels.

Below are example atmospheric parameters that produce the problem.

Using robust statistics on the output photons indicates that most photons land in a reasonable location, with a set of outliers

I can also provide a full example imsim config that produces this set of atm parameters

{'L0': [17.174188473465833,
        17.174188473465833,
        17.174188473465833,
        17.174188473465833,
        17.174188473465833,
        17.174188473465833],
 'altitude': [1.4302596307983027,
              5.348763607590827,
              11.208179050816556,
              18.24311343973453,
              24.467319796805825,
              30.175959028659584],
 'direction': [coord.Angle(5.023031335581778, coord.radians),
               coord.Angle(5.4138391122102645, coord.radians),
               coord.Angle(5.213667959962739, coord.radians),
               coord.Angle(5.043782770080258, coord.radians),
               coord.Angle(5.363918865113014, coord.radians),
               coord.Angle(5.075185823976301, coord.radians)],
 'r0_500': 0.08161434310901314,
 'r0_weights': [5.531557367962777e-13,
                4.3652382380659045e-14,
                1.326839454952542e-13,
                4.126186570683441e-14,
                2.344518516479255e-14,
                2.5838659853736628e-15],
 'rng': galsim.BaseDeviate('1778941255 1431167650 1282975060 ... 4184287197 2656677445 2301135708'),
 'screen_scale': 0.1,
 'screen_size': 800.0,
 'speed': array([10.94986395, 17.84175431, 37.29034581, 42.9768206 , 19.97968586,
       10.84201688])}

[rmjarvis]

I think I need more details to reproduce this. I tried to turn this set of kwargs into a unit test, but it doesn't fail yet. Can you try to adjust this based on other details of your run until it does fail? Then Josh or I can have a better shot at diagnosing and fixing it.

My first attempt is on branch '#1280' here:
https://github.com/GalSim-developers/GalSim/blob/%231280/tests/test_phase_psf.py#L1557

[rmjarvis]

Note: If the issue is just the rng, the repr isn't sufficient to serialize it, since the repr string is really long, and that usually gets annoying to see long rngs all over the place. You can get a full evalable string for the rng seed with rng.serialize().

[esheldon]

I can't get this to fail even either, even with the full serialized rng.

It's complex because the same atmosphere gets used for all the objects in the image, and there are many different star SEDs being used. I see the problems for some fraction of the objects.

I could give you something that will reproduce using an imsim run. I would give you the input data. You would also need to use my hacked up imsim so it will printout offset statistics as well as when it hits a nan.

[rmjarvis]

Maybe you could print out the repr of the object whose drawImage call produces the nans? And also any other kwargs being sent to drawImage. That should have enough information to reproduce it I think.

[esheldon]

I'll try that. I'll also need to repr the photon_ops in that case, since that is how the PSF is implemented

[esheldon]

It is not easy to get a perfect reproduction. There must be some internal settings from the imsim run that affect the results.

Also note I had to turn off sensor, which made the run too slow. It is not that slow in the imsim run so, again, there must be some settings I'm not aware of.

Some objects from imsim do not have a full repr and cannot be pickled, so I tried to reconstruct things by hand when necessary.

I've attached a tarball that has data as well as a file test.py that can be run with pytest -vv test.py and will see failures of the tests for nan and large offsets.

I wasn't sure what an acceptable range for offsets should be. In this example the max offset I saw was 600_000 pixels (120_000 arcseconds). I have a check right now at 10_000 pixels, but we should adjust it.

nan-example.tar.gz

[esheldon]

I'm aware that this can't be incorporated into the galsim test suite due to use of imsim, but it wasn't obvious to me how to reconstruct the photon_ops without using imsim code

[rmjarvis]

Thanks. I'll work on minimizing the example to make a unit test from it.

[rmjarvis]

So the first thing I'm noticing here is that you have RubinDiffractionOptics() in the photon_ops list. So that explains why you are getting large offsets. You should expect large offsets from the diffraction spikes.
And I suspect that the NaNs are also due to that. If I remove that item from the photon_ops, the NaNs go away too. I'll try to track down where the NaNs first show up. But I think this is most likely an imsim bug, not a GalSim one. (I won't close this until I confirm though.)

[rmjarvis]

@jmeyers314 Josh, I tracked this down to a batoid function. Specifically, refract in trace.py. But it happens in the C layer, so I'm passing the baton to you now.
One of the ray vectors has the following values before the _batoid.refract call:

x,y,z,vx,vy,vz,t,wave,flux,vig,fail = 1.4136048029448713 -1.26631532296121 0.5702413527436976 -0.1632631352399052 0.34430459097839733 0.5688448339442601 17.841282831343328 5.356478398425508e-07 0.3783200003476515 True False

After refract:

x,y,z,vx,vy,vz,t,wave,flux,vig,fail = 1.619778407831917 -1.701113531367937 -0.20811303793564317 nan nan nan 16.578452709992995 5.356478398425508e-07 0.3783200003476515 True False

More details:

surface =  Sphere(-13.36)
ct =  CoordTransform(CoordSys(array([0.        , 0.        , 4.34206865]), array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]])), CoordSys(array([0.        , 0.        , 4.40206865]), array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]])))
m1 =  SellmeierMedium((0.69618302, 0.407925877, 0.897464057, 0.00467926519, 0.0135122244, 97.9323636))
m2 =  ConstMedium(1.0)
rv =  RayVector(array([0.06317926, 0.05853737, 0.0911232 , ..., 0.0587624 , 0.05450346,
       0.09584311]), array([-0.03250948, -0.03895119, -0.0603161 , ..., -0.03323284,
       -0.03971611, -0.03526725]), array([0.00079655, 0.00078003, 0.00188413, ..., 0.00071907, 0.00071758,
       0.0016456 ]), array([ 0.05435319,  0.08560945, -0.12966972, ...,  0.08319322,
        0.1107581 , -0.15863831]), array([-0.10171135, -0.0591399 ,  0.08145143, ..., -0.09638325,
       -0.05337458, -0.08124123]), array([0.67436141, 0.67677561, 0.66551091, ..., 0.67180015, 0.66982595,
       0.66125685]), array([17.64212032, 17.64272894, 17.59542174, ..., 17.64550459,
       17.64643467, 17.60539304]), array([5.12374699e-07, 5.37952130e-07, 4.67085877e-07, ...,
       4.96827704e-07, 4.16832107e-07, 5.44127040e-07]), array([0.37832, 0.37832, 0.37832, ..., 0.37832, 0.37832, 0.37832]), array([False, False, False, ..., False, False, False]), array([False, False, False, ..., False, False, False]), CoordSys(array([0.        , 0.        , 4.34206865]), array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]])))

So it looks like the diffraction sent one ray out past the edge of the focal plane. It has vv large x,y,z values and it's marked as vignetted. That should be fine, but something in here is turning vx,vy,vz into nans. Hopefully that will be enough for you to see what might be wrong.

[jmeyers314]

Thanks, I'll take a look. Note though that if it's marked vignetted, it shouldn't end up in the final image (flux should be set to 0 inside ImSim somewhere).

[rmjarvis]

Probably so. I think Erin was doing an empirical centroid with them just took the mean of the x and y values (which had gotten turned into NaNs subsequent to this step). So Erin, probably if you exclude the flux=0 photons, the centroids will be finite. But they still aren't doing what you want, since the centroid of the diffraction spikes are going to be quite noisy.

[jmeyers314]

Looks like this particular refraction is actually a case of total-internal-reflection, which I haven't properly addressed in batoid. In lieu of actually reversing the order of the trace, I'll add code to batoid to mark such cases as failed for now.

[esheldon]

All the im.photons.flux are set to nan in all cases