Simplify draw procedure for very faint objects

[rmjarvis]

This PR implements a proposal I made here to simplify the draw procedure for very faint things.

For the objects that are only going to have of order a 100 photons or less, we mostly just care that the flux is there. Not that they have a realistic wavelength spectrum and careful treatment of sensor effects that we care about for brighter objects that will actually be used for science.

There are two tunable parameters that we can play with to trade fidelity vs speed:

1. max_flux_simple is the maximum flux (in any band) we allow for objects that will get these simplifications. Default = 100
2. sensor_limit is the minimum fluctuation already existing in the image that would force the use of the full SiliconSensor model. Default = 200

Only if the flux is less than max_flux_simple and there aren't any pixels yet above sensor_limit do we turn off the sensor effects.

But the SED is simplified either way if the flux is less than max_flux_simple, since that only impacts very subtle things like the exact shape of the DCR effect and the depth distribution of the converted electrons in the silicon.

[cwwalter]

Just confirming my understanding:

Surface operations work even without the silicon model right, so the effect for DCR is just that the SED is wrong but it will still smear the object out in someway.

[rmjarvis]

Surface operations work even without the silicon model right, so the effect for DCR i

Yes. They require photon shooting, but they are independent from the type of sensor used.

[cwwalter]

Wait... I'm not sure we want to put this into master yet...

Sorry, it's late here but I still have questions about a particular piece of it.

[cwwalter]

Need to go to sleep here right now but one question I wonder if it is easy for you to answer by doing a run Jim:

If you change the new parameters (max_flux_simple, and sensor_limit) so that all of the changes except the skipping of the sensor are still applied( i.e. trivial SED, new center, new RA etc etc) how much of a speed up do we have?

I mean, I am interested in knowing if we have the same behavior as before for always using the silicon, but with the other speedups applied, how much of the ~x2 gain do we lose?

[cwwalter]

And a question for Mike: we are also talking about making the carpet of dim objects by drawing them not as individual objects but drawn from a pre-realized sky addition.

Would you see this as being a class of objects between the fully simulated objects and the ones in the carpet? Or do you think this change makes things fast enough it would negate the need for the carpet technique?

[rmjarvis]

I think you can almost consider the simple sed and no sensor as one way to implement the carpet idea. I think these really faint things aren't now much slower than what you would probably be able to do with the carpet idea. Except here we get the variable PSF right. (The carpet would probably require approximating the PSF as constant across a CCD.)

[rmjarvis]

For the always-use-the-sensor test, you can set sensor_limit=0.

I predict it will lose quite a lot of the gain. From the latest cProfile output, if I've run the numbers right, then bumping the fraction of objects that have to run the accumulate function will slow things down by about 25% overall. Said differently, if it was a factor of 2 (50%) faster, then now it will only be 37% faster.

[cwwalter]

I think you can almost consider the simple sed and no sensor as one way to implement the carpet idea. I think these really faint things aren't now much slower than what you would probably be able to do with the carpet idea. Except here we get the variable PSF right. (The carpet would probably require approximating the PSF as constant across a CCD.)

Yes, that it what I was wondering. This would be really good if we could do it instead. I think dealing with each object separately has many advantages that we would lose in the other scheme.

OTOH, there is an issue that I am concerned about that I would like to understand a bit better. From my previous experience changes like this can cause issues in the sense that we may wonder if subtle things like PSF errors have been introduced and (even if it is not the case) when trying to understand something that issue or bug we will notice, we will wonder if this is what is causing it.

With something like the carpet scheme it has the advantage that we need to make an external object drawing decision based on an a criteria set on the input quantities (like say magnitude). So, this has the advantage that, after matching, we can exclude or separately examine that sub sample of objects. With a procedure like this where we do a test based on the running condition we can't, and I think that is problematic.

So, I was thinking about approaches to this issue since, if we could either replace the carpet idea with this (which has a several storage, psf and tracking issues), or if it turns out (as is probable) that it is a really a big part of the speed up we will probably want to do it. One approach we could take is to add some more information to the centroid file to classify how objects were handled. I think there are at least 4 categories now (which is one of the things that is worrying me):

• Fully simulated
• Simple/Faint
• Skipped
• (in future maybe) carpet.

Adding an entry which tracked that would at least allow us to isolate sub-samples for study purposes. I was already worried that there was no record in the files for the skipped cases.

[cwwalter]

Two other tangential issues from reading the code and thinking about this:

• Mike mentioned sticking lots of things/information in HDUs in FITS files in the meeting yesterday. When I started thinking about extending the centroid files (and also remembered some of the atomic write issues we had on Theta) I wondered if we could save both disk space (since we could use binary info) and number of files by putting the information in a HDU which would travel with the FITS file instead.

This has a few disadvantages though: the FITS files would be bigger, you would need a standalone program to either read them or produce text files, and if you couldn't just store the centroid files separately for analysis purposes (unless you created them separately after the fact).

• One of the issues in the sims_galSimInterface code is that the spit between the sensor model and non-sensor model case is awkward. You have to put some of the same things/almost duplicate code for both cases. That is why there is currently not yet DCR in the non-sensor model case now for example. But, this PR from Mike is now making decisions in the sensor version of the code about when to skip actually using the sensor model.

Currently they are in different classes, and I was too tired last night to go through carefully but I wonder with this framework in place if we could now combine the cases and solve some of these issues (like having to turn the DCR on in both places).

[jchiang87]

I would note that since the centroid files contain the model fluxes and not the realized fluxes, they can really be generated at any time, i.e., they don't have to be created at the same time that the full simulations are running. This would give us a lot more flexibility in how they are handled. In fact, it would be useful to compute them before the full simulation since that info could be used to predict the runtime of the full simulation.

[cwwalter]

I would note that since the centroid files contain the model fluxes and not the realized fluxes, they can really be generated at any time, i.e., they don't have to be created at the same time that the full simulations are running. This would give us a lot more flexibility in how they are handled. In fact, it would be useful to compute them before the full simulation since that info could be used to predict the runtime of the full simulation.

Interesting/good idea...

I think we would want to be able to do it from the instance catalogs (as opposed to the process that made them) so that people could make them for any instance catalog that they made. Perhaps it could be part of the suite associated with the binary instance catalogs that can also do a translation back and forth to text.

But could we really know ahead of time if (for example) we had decided to treat a galaxy as faint? That will actually depend on the running conditions right? I would like to be able to mark that.

I had been thinking about adding some actual realized quantities to the file, but I suppose that could go elsewhere if we did this.

[jchiang87]

@danielsf If you haven't done so already, could you please pass this through Jenkins? I'm doing all my work from docker images which have been stripped of some test data in lsst_distrib so I can't run all of the tests myself.

[danielsf]

Starting the job now

[danielsf]

@jchiang87 This branch passed Jenkins

[danielsf]

I'll merge it at the end of the day today, unless someone objects

[cwwalter]

I think this code currently looks good. And it should be matched with Jim's latest: https://github.com/LSSTDESC/imSim/pull/181

But, as I said above, I don't think I am comfortable with using the sensor_limit etc parameters in non-zero mode yet in non-test running. This actually corresponds to some variable magnitude cut-off to do full object handling. Especially if this replaces the "carpet" (which would be nice) I think we need to have some discussion about where that cut should be and also do some tests to understand how the speed up scales with magnitude cut and what fraction of the speed up is do to this vs the other fixes.

Then, since this is something that can only be determined at run time (unlike a predefined magnitude cut) and it will depend on the sky level and nearby objects, I think we should mark these galaxies (I think the centroid file is really the only place we can do it?) so we can know after the fact which ones were handled this way.

Also, I wonder how much extra time we are now saving completely skipping the "zero" flux cases if we are doing this too.

Finally if I do make a handmade small instance catalog with dim sources and I want to see their effect; I might not want this turned on since I may not care about the run time.

So, I think I am proposing merging this but with the default set in line 244 for not applying the cut with the ability to turn it on in a patch like Jim's.

[rmjarvis]

Especially if this replaces the "carpet" (which would be nice) I think we need to have some discussion about where that cut should be and also do some tests to understand how the speed up scales with magnitude cut and what fraction of the speed up is do to this vs the other fixes.

To be clear, the "carpet" mode is equivalent to sensor_limit = inf. So this is already significantly more conservative than that. The carpet mode would basically never use the sensor model properly. We can move to that if you prefer to have all faint objects treated equivalently, which would speed things up (probably only slightly) more.

[rmjarvis]

Also, I wonder how much extra time we are now saving completely skipping the "zero" flux cases if we are doing this too.

That's not a huge effect. In the latest profile that Jim sent out, there were 196399 objects requested to be drawn (drawObject), and 160871 made it to actually drawing something (drawImage). So about 18% of objects apparently had 0 flux and were skipped.

[rmjarvis]

I am proposing merging this but with the default set in line 244 for not applying the cut with the ability to turn it on in a patch like Jim's.

Done. e34bed7

[danielsf]

@cwwalter I will let you make the final call when this gets merged.

[cwwalter]

Done. [e34bed7]

Thanks! Could we also change the max_flux_simple for the default? I realize this will work as is as long as there is some sky background but there is a case I can think of where we will want both: doing detailed studies of the sensor model. So, for example, in the past in PhoSim I made spots as a function of intensity and position in the sensor model to track the BF effect with all other effects turned off (including the atmosphere). I think we want people to be able to do that unless we explicitly turn them on in the config file.

@danielsf after that change I think we should go ahead and merge this.

[cwwalter]

To be clear, the "carpet" mode is equivalent to sensor_limit = inf. So this is already significantly more conservative than that. The carpet mode would basically never use the sensor model properly. We can move to that if you prefer to have all faint objects treated equivalently, which would speed things up (probably only slightly) more.

No, I agree with you that this this is better. No external storage, and proper PSF handling.

My ideal case is that we only do this (no carpet, no object skipping) and that we have a way of marking them so that we know the cases/magnitude for when this behavior kicks in.

I guess we will get some more speed up when with the new binary instance catalog if we can group all galaxy components together so we don't have to do the sizing operations multiple times.

I think timing a run with our guess for the correct parameters and without object skipping would be good. It would tell us if the overall performance is OK with only your new code.

[rmjarvis]

Could we also change the max_flux_simple for the default?

To what? I'm confused. The default now is just to use a simpler SED for the faint things, not turn off the sensor effects. So the rest of your comment doesn't really make sense in relation to changing max_flux_simple.

I added a bit more to the doc string to hopefully make this clearer, but please let me know if you nonetheless still want the default for max_flux_simple to be changed (and to what?).

[cwwalter]

Hi Mike, I'm sorry I didn't manage to get to this yesterday.

OK I think there are two issues: For:

The default now is just to use a simpler SED for the faint things, not turn off the sensor effects. So the rest of your comment doesn't really make sense in relation to changing max_flux_simple.

Going back and looking more carefully; I think I was confused by the comment:

                    # For faint things, only use the silicon sensor if there is already
                    # some significant flux on the image near the object.
                    # Brighter-fatter doesn't start having any measurable effect until at least
                    # around 1000 e-/pixel. So a limit of 200 is conservative by a factor of 5.
                    # Do the calculation relative to the median, since a perfectly flat sky level
                    # will not have any B/F effect.  (But noise fluctuations due to the sky will
                    # be properly included here if the sky is drawn first.)

Since this in fact says that for faint things the sensor will be turned off unless there is something bright around it.

But, in fact it depends on the value of sensor_limit (which is just referred to as a constant of 200 above) and the code

            if np.max(image.array) > np.median(image.array) + sensor_limit:
                        sensor = self.sensor[detector.name]

does what you say since if sensor_limit is zero this is always true and the sensor will be used. But I was thinking that the sensor could still be turned off (which is not correct).

For the 2nd thing:

More generally, I was arguing that I preferred the default to be that every object be completely simulated including a full SED. I think is is fine to have the default config file in imSim call this with these parameters set to speed things up, but this is code most people would never look at (since it isn't even in the imSim code base) and it will do something drawing these objects that is not obvious to non-experts of the code. I can't exactly anticipate the test people will do that would bite them but I can imagine doing something at low flux level and then adding things together and not seeing the SED you put in.

In a DC2 type context, I think we certainly want this on, but for stand alone smaller applications we might not and I would rather people explicitly understand they are turning on this speed up and what the consequence will be. So, to this we would make max_flux_simple=0 too I think, correct?

One documentation comment I noticed when going through this again:

Since the non sensor based class comes first the comments:

        @param [in] max_flux_simple is ignored here.
        @param [in] sensor_limit is ignored here.

Are a bit confusing since, if you are reading linearly, you don't know what they are. Maybe adding: 'They are used in the GalSimSiliconInterpeter drawImage.' would be nice.

Thanks.

[rmjarvis]

OK. Done.

[cwwalter]

OK thanks. @danielsf I think this is ready to merge whenever you are ready.