Skip to content

Refactor image/stamp workflow#424

Merged
welucas2 merged 63 commits intomainfrom
u/g-braeunlich/refactor-image
Jun 5, 2025
Merged

Refactor image/stamp workflow#424
welucas2 merged 63 commits intomainfrom
u/g-braeunlich/refactor-image

Conversation

@g-braeunlich
Copy link
Contributor

No description provided.

@g-braeunlich g-braeunlich self-assigned this Sep 28, 2023
@g-braeunlich g-braeunlich force-pushed the u/g-braeunlich/refactor-image branch 2 times, most recently from e58012d to c145f0c Compare October 10, 2023 12:34
@g-braeunlich g-braeunlich force-pushed the u/g-braeunlich/refactor-image branch from c145f0c to 1d7fefc Compare October 30, 2023 09:16
@g-braeunlich g-braeunlich force-pushed the u/g-braeunlich/refactor-image branch 2 times, most recently from acbac83 to 9f49df7 Compare November 30, 2023 12:45
@welucas2
Copy link
Collaborator

The latest commit aligns photon objects in the photon pooling pipeline with their positions in the original pipeline.

@welucas2 welucas2 self-assigned this Apr 29, 2024
@welucas2
Copy link
Collaborator

Almost there, but blocked by GalSim-developers/GalSim#1284.

@welucas2 welucas2 force-pushed the u/g-braeunlich/refactor-image branch from daae77c to 7290a11 Compare May 24, 2024 12:11
@welucas2 welucas2 force-pushed the u/g-braeunlich/refactor-image branch from bdaadae to b786cb0 Compare June 5, 2024 13:10
@welucas2 welucas2 marked this pull request as ready for review June 5, 2024 13:32
@welucas2
Copy link
Collaborator

welucas2 commented Jun 5, 2024

This is finally open for review -- all comments welcome.

Copy link
Contributor

@rmjarvis rmjarvis left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks William. I have a bunch of comments, but this is definitely a good start. If it would be helpful to talk through any of this on a zoom call, let me know.

# If using photon pooling, FFT and photon objects are batched separately.
# There will probably be few enough FFT objects to do in a single batch.
# Bright photon objects are treated in all photon batches, but at 1/nbatch_photon of
# their flux. Faint photonobjects are placed at their full flux in random batches.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
# their flux. Faint photonobjects are placed at their full flux in random batches.
# their flux. Faint photon objects are placed at their full flux in random batches.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for spotting this. Fixing in the next commit.

# These parameters only affect LSST_PhotonPoolingImage images.
# The batch numbers can be changed if desired.
nbatch_fft: 1
nbatch_photon: 10
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this should be the same parameter as the simple nbatch. We can just document that with photon pooling, nbatch refers only to the photon-shooting objects. Then nbatch_fft can be an optional parameter that defaults to 1. (Maybe not even worth mentioning here if most users wouldn't change it?)

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree with this. Nice to keep the number of parameters down. nbatch_photon will be no more, and the photon pooling will be controlled directly with nbatch.


class LSST_ImageBuilder(LSST_ImageBuilderBase):
"""This is mostly the same as the GalSim "Scattered" image type.
So far the only change is in the sky background image."""
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This doc is probably obsolete. There are quite a few differences now.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Happy to change this. Can you suggest something more suitable?

'x' : { 'type' : 'Random' , 'min' : xmin , 'max' : xmax },
'y' : { 'type' : 'Random' , 'min' : ymin , 'max' : ymax }
}
set_config_image_pos(config, base)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Rather than making this a free function in a different file, this should probably be a private method of the base class, which both subclasses can call.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed, doing this.

base["stamp"]["photon_ops"] = []
photon_ops = base["stamp"]["photon_ops"]
shift_op = {'type': 'Shift'}
shift_index = next((index for (index, d) in enumerate(photon_ops) if d["type"] == "Shift"), None)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't love this implementation. I get why something like this is now required. But I think I'd rather have RubinOptics and RubinDiffractionOptics just have an optional parameter to apply this shift or not. Then when running this class, base[image_pos] is defined and they can apply the shift based on that.

In the photon pooling class, once we've pooled all the photons, we should set the object-specific values (image_pos, sky_pos, maybe others) that apply to a particular object to None. Then the RubinOptics and RubinDiffractionOptics photon ops would see that the image_pos isn't defined, so not apply any shift.

This feels to me more elegant and possibly more efficient.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should now be done -- see RubinOptics and RubinDIffractionOptics. If an optional shift_photons = True is passed, the shift is performed.

},
"psf": {"type": "Convolve", "items": [{"type": "Gaussian", "fwhm": 0.3}]},
"stamp": {
"type": "LSST_Silicon",
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should be the new stamp type when image_type is LSST_PhotonPoolingImage

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, photon pooling now needs to use the new stamp type rather messing with base. The test sets the correct stamp type.


def test_lsst_image_photon_pooling_pipeline():
"""Check that LSST_PhotonPoolingImage batches objects as expected and renders objects at the correct positions."""
run_lsst_image("LSST_PhotonPoolingImage")
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I feel like we probably want a few additional tests in this pipeline. E.g. that the FFT and faint objects only got drawn once. And the others got drawn nbatch times.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've written a set of tests for the partitioning and batching process. As mentioned above, it also confirms that total flux is conserved for photon objects,

help="Similar to -k of pytest / unittest: restrict tests to tests starting with the specified prefix.",
default="test_",
)
args = parser.parse_args()
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

FWIW, when working on stuff, I usually just put the name of the test function I'm working on at the start of this name == '__main__' block, followed by quit(). Then remove it before committing. Seems easier than rolling all of this argparse stuff, but 🤷.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is also done in test_diffraction_fft.py, but I'm happy to change it if you'd like me to.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess if you find it useful, go ahead and keep it. But if you were just copying from that file, I don't think this is necessary to keep.

expected_x_pic_center = 564.5
expected_y_pic_center = -1431.4
expected_x_pic_center = -989.5971378245167
expected_y_pic_center = -3840.3512012842157
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why are these (and subsequent tests) changing so much? This isn't close.

Regression tests are supposed to make sure we don't change the functionality of existing code. I know the implementation of these photon ops changed, but I think we want the tests to remain the same or at least very very close. What happened here?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is still a big issue. I'm 99% convinced the error is in the original test itself. My current best guess is that with the old code it may not have been internally self-consistent around the WCS used, while the new XyToV should enforce self-consistency. But it still needs to be looked at.

"bandpass": galsim.Bandpass('LSST_r.dat', wave_type='nm'),
"wcs": galsim.PixelScale(0.2),
"wcs": wcs,
"current_image": galsim.Image(1024, 1024, wcs=wcs)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why was this change required?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The change to the wcs field isn't actually necessary, it looks. The addition of current_image is needed for deserialize_rubin_optics, deserialize_rubin_diffraction_optics and deserialize_rubin_diffraction which uses it to get img_wcs, required by the new version of XyToV.

@welucas2
Copy link
Collaborator

welucas2 commented Nov 7, 2024

The latest set of commits should address most of the comments above. I'll respond to those individually. A couple do still need to be looked at.

@rmjarvis
Copy link
Contributor

Are you waiting on me for anything here? I thought you still had some more comments to address, so I was waiting for that to be done. But if it would be helpful for me to review this again at this point, let me know.

@welucas2
Copy link
Collaborator

Not just yet -- I have a couple more commits on the way with fixes and tests for the batching. I'll push these later today or tomorrow, and then I think that should be everything ready for review again.

@welucas2 welucas2 force-pushed the u/g-braeunlich/refactor-image branch from 46461fb to eeb67dc Compare November 20, 2024 16:54
@welucas2
Copy link
Collaborator

I think this is ready for re-review now, though CI is failing with an error in test_stamp_bandpass_airmass from numpy coming through GalSim. I saw the same error a couple of weeks ago on an ARM system when using GalSim 2.6 which introduces its own utilities.least_squares(). I wasn't able to investigate further before the system went down and only came back up yesterday.

@welucas2
Copy link
Collaborator

I think there is actually one more commit incoming. I've realised that in edge cases where nominally bright photon objects' phot_flux falls below nbatch, we can end up with batches containing objects with 0 flux. I think I have a fix, by having these objects be treated as faint objects for batching purposes. I'll test this and if all is well push it tomorrow.

@welucas2
Copy link
Collaborator

OK, that should be it -- code open for review!

@welucas2 welucas2 requested a review from rmjarvis November 22, 2024 17:05
Copy link
Contributor

@rmjarvis rmjarvis left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks William. I think we're close.

I'd still like to have a separate PR that just changes the tests of the expected centers. I think the way to do that is to change the test to mostly your new code, but when calling (the old) XyToV, use sky_pos = img_wcs.toWorld(image_pos) and local_wcs = img_wcs.local(image_pos). I'm not sure what image_pos is exactly, but hopefully you can follow Geri's code to see what the implicit image_pos is there.

offset_adjustment: a PositionD by which to offset each object
"""
for stamp in stamps:
offset = offset_adjustment
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Probably should just call this parameter offset from the start.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks! Doing this now.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And now obsolete with the big simplification in offsetting the photons by stamp_center. I've removed the two static methods in photon_pooling.py dealing with the offset and it's now dealt with in stamp.py just using stamp_center.

Returns:
images: Tuple of the output stamps. In normal PhotonPooling usage these will actually be
PhotonArrays to be processed into images after they have been pooled following this step.
current_vars: Tuple of variables for each stamp (noise etc).
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These should be the current noise variance in each image. var=variance, not variable.
It's not really relevant for most imsim usage, but if you draw a RealGalaxy, and especially if you then whiten it, the image gets some noise from that, so this mechanism helps GalSim from adding more noise than requested when adding the rest of the image noise.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for pointing this out. Is the actual usage of current_vars ok?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes. Just the doc string is wrong.

Comment on lines +313 to +315
batches = [
[dataclasses.replace(obj, phot_flux=np.floor(obj.phot_flux / nbatch)) for obj in phot_objects]
for _ in range(nbatch)]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A nifty (simpler?) way to get the total flux to come out right:

Suggested change
batches = [
[dataclasses.replace(obj, phot_flux=np.floor(obj.phot_flux / nbatch)) for obj in phot_objects]
for _ in range(nbatch)]
batches = [
[dataclasses.replace(obj,
phot_flux=(obj.phot_flux*(i+1))//nbatch - (obj.phot_flux*i)//nbatch)
) for obj in phot_objects]
for i in range(nbatch)]

Then you don't need the block below this.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Very nice -- thanks! Putting this in now.

imsim/stamp.py Outdated
def build_obj(stamp_config, base, logger):
"""Precompute all data needed to determine the rendering mode of an
object."""
stamp_type = stamp_config.get('type','LSST_Photons')
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think we want to have a default here. Rather we should make sure the user has correctly specified this type. So if stamp_config.get('type', None) != 'LSST_Photons':... then raise some kind of explanatory errror.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks -- doing this now.

imsim/stamp.py Outdated

max_flux_simple = config.get('max_flux_simple', 100)
faint = self.nominal_flux < max_flux_simple
bandpass = base['bandpass']
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

duplicate

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

D'oh! Thanks for spotting -- fixed.

imsim/stamp.py Outdated
img_pos = base["image_pos"]
obj_offset = base["stamp_offset"]
image.photons.x = image.photons.x + img_pos.x - obj_offset.x
image.photons.y = image.photons.y + img_pos.y - obj_offset.y
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think all this offset stuff could be simpler. Here you add img_pois - stamp_offset. Then later you add 0.5 for even sized images. So the net adjustment is

img_pos - stamp_offset + (0.5 if size%2==0)

I'm pretty sure this combination is equivalent to base["stamp_center"] based on the calculation of stamp_offset here. So I think you can just add stamp_center here and skip the later adjustments. Which honestly makes sense. These photon arrays are meant to be relative to the image center, which here is just the stamp center.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Very nice call. I've just tested and this gives identical output. Making this change now.

# positions the photons before they are pooled to be processed together.
if self.shift_photons and self.image_pos is not None:
photon_array.x += self.image_pos.x
photon_array.y += self.image_pos.y
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Per my other comment, these should presumably also shift by stamp_center, not image_pos. Likely not a big difference, but might as well use the correct shift.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this also be done at the end of applyTo? e.g.

imSim/imsim/photon_ops.py

Lines 111 to 112 in a4c2d97

photon_array.x -= self.image_pos.x
photon_array.y -= self.image_pos.y
in main as it is now.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes. I think so.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is the last thing on my list. It means changing the API to RubinOptics and RubinDiffractionOptics. I think normal use is always through photon op registrations via the photon_op_type decorator and deserialize_* functions, which wouldn't need to change, but at least the tests in test_photon_ops.py do directly init the photon operators themselves. Is that ok to change?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think if we preserve the user-facing code (i.e. how the config file specifies things), then it's ok to make changes to the implementation, even if these are technically API changes. But perhaps @jchiang87 or @cwwalter might want to weigh in about how strictly we want to preserve backwards-compatibility of the Python classes/functions.

help="Similar to -k of pytest / unittest: restrict tests to tests starting with the specified prefix.",
default="test_",
)
args = parser.parse_args()
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess if you find it useful, go ahead and keep it. But if you were just copying from that file, I don't think this is necessary to keep.

"""

def __init__(self, local_wcs, icrf_to_field, sky_pos):
def __init__(self, local_wcs, icrf_to_field, img_wcs):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this doesn't need local_wcs anymore. If we're breaking the API (with the change in the last parameter), might as well clean it up and remove local_wcs as well.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done. This also leaves some of the photon_velocity methods in RubinOptics, RubinDiffraction and RubinDiffractionOptics with unnecessary local_wcs and, now I look at them, also rng. Would you prefer those left alone or removed as well?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Might as well clean up all the obsolete API. If any of them seem like user-facing code, we should make this a deprecation rather than completely break old code. But I think most of these functions are internal, not user facing, so it should be ok.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I forgot to mention this in the meeting! Yes, everything looks to be internal to the photon ops themselves. rng can't be removed though, as it turns out, or it breaks RubinDiffractionOptics. Otherwise, tidied up as requested.

@welucas2
Copy link
Collaborator

welucas2 commented Dec 6, 2024

Assuming tests pass, that should be everything above. I've included changes to RubinOptics to use stamp_center rather than image_pos, changing the tests to fit, but let me know if you don't like that and I'll revert it.

I'll move on to seeing if I can make the old photon ops tests match the new ones as suggested.

@welucas2
Copy link
Collaborator

I've just created PR #484 which corrects the photon op tests to be self-consistent. There were actually a couple of inconsistencies in the tests here too, but working through them now means everything should be in alignment pre- and post-photon pooling.

I hope this means that once #484 is merged in, this can be, too.

Copy link
Collaborator

@jchiang87 jchiang87 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi William,

In comparing how this branch works wrt the code in main, I saw a few items that may need to be addressed.

# If using LSST_Image (non-photon pooling), nbatch is the overal number of batches in
# which to draw the objects.
# The default number of batches is 100, but you can change it if desired.
nbatch: 100
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be good to clarify that nbatch is also used by the photon pooling code, not just by LSST_Image.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks -- doing this now.

imview._shift(-center) # equiv. to setCenter(), but faster
imview.wcs = PixelScale(1.0)
if imview.dtype in (np.float32, np.float64):
sensor.accumulate(photons, imview, imview.center)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there ever a situation where we would want to set resume=True here, e.g., similar to this case in galsim?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure -- it looks like it can boost performance somewhat when using a SiliconSensor by skipping an initial calculation in the accumulate method. I don't think there's any disadvantage to using resume=True, but it might be worthwhile discussing this further in tomorrow's call to be certain.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Making this change now.

Copy link
Collaborator

@welucas2 welucas2 Jan 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, though the tests pass it looks like my implementation has problems in practice as it's a different image view being used each time. I'm trying a fix now by creating the image view outside the loop over photon pools and reusing it.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should be fixed with the latest commit.

imsim/stamp.py Outdated
image=image,
wcs=base['wcs'],
sensor=NullSensor(), # Prevent premature photon accumulation
photon_ops=psfs,
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this be

photon_ops=photon_ops,

?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You're right, it should be photon_ops, but the code above which creates the initial photon_ops list from the config is incorrect -- thanks for noticing something's not right here. It should be contain the PSF and the BandpassRatio photon op (only if the SED needs to be reweighted), which need to be dealt with here. All the other photon ops are dealt with later in the pipeline and applied to the whole photon batch. I'm fixing this now.

save_photons=True)
stamp_center = base['stamp_center']
image.photons.x = image.photons.x + stamp_center.x
image.photons.y = image.photons.y + stamp_center.y
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The non-photon pooling code sets

base['realized_flux'] = image.added_flux

at the corresponding point in the code so that those values appear in the output.truth files. What should it be set to here so that the equivalent info is written out?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a tricky one. Per object quantities are diluted with the use of photon pools. If we did the same thing here in the stamp, we'd firstly only be dealing with an object flux which is a fraction ~1/nbatch of the total flux, complicated by different pools having different integer values of the flux in order to conserve the total. Furthermore it might have been possible to track and sum the realised flux per object across the pools, but the accumulation now happens at the pool level. We can tell how much of the pool's flux is added to the image, but the information about which object contributed which photon is not retained at that point. I'm not sure if someone else has a suggestion or is aware of something that would help. I'll have a think about this but I think we might need to chat about this in tomorrow's call.

@welucas2 welucas2 force-pushed the u/g-braeunlich/refactor-image branch from bf6cb43 to 72cc687 Compare January 17, 2025 11:54
@welucas2 welucas2 force-pushed the u/g-braeunlich/refactor-image branch from f597359 to 79d7eac Compare May 9, 2025 11:50
@welucas2
Copy link
Collaborator

welucas2 commented May 9, 2025

Rebased onto #485's fix to the CI, tests are now passing. @rmjarvis, @jchiang87, are you happy to merge?

# nbatch is used for LSST_Image and LSST_PhotonPooling runs. In the former it determines
# the number of batches into which to all objects are placed, while in the latter it is
# used only with photon objects (bright and faint).
self.nbatch = params.get('nbatch', 100)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If we're doing subbatches, I think 100 is probably overkill for the main batch number. 10 or 20 should be sufficient for the science needs of getting BFE correct. Maybe adjust the distribution of these two. (E.g. default nbatch=20, nsubbatch=50 to keep 1000 total.)

Also (related) nsubbatch isn't included in the config file, which kind of serves as our documentation of this stuff. Probably worth adding there.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What about for LSST_Image runs -- are we happy changing the default for people who will continue using that process? I could have it use a different value depending on image type, 100 for LSST_Image, 10 or 20 for LSST_PhotonPoolingImage.

I set nsubbatch in the new photon pooling config (imsim-config-photon-pooling.yaml) since it isn't used for LSST_Image runs. I can move it to imsim-config.yaml if you'd like, or document it where it is now.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Our template config file for LSST_Image sets nbatch explcitly, so maybe leave it 100 there, but change the comment to say that the default is 20, and that there is a speed/memory tradeoff involved, so when memory is at a premium, increasing it to something like 100 is generally warranted.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks. I've switched things around to move nsubbatch from imsim-config-photon-pooling.yaml into imsim-config.yaml, and added some docs, and used your suggested defaults of nbatch=20 and nsubbatch=50 for photon pooling (and keeping nbatch=100 for LSST_Image). I might push a little for a default of nbatch=10 for photon pooling with which runs take about 60% the time of nbatch=20, but at least it's something that can be easily tuned.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm fine with a default of 10, rather than 20. I think almost all use cases will not care about the difference between these. 20 is more conservative, but I doubt 10 will ever be insufficient.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've changed again to 10 and noted this in the config comment.

op.applyTo(photons, local_wcs, rng)
# Shift photon positions to be relative to full_image.center.
# This is necessary as all photons will now be in the pool together, and there
# is no longer any way to distinguish which belong to which object.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's more that this is necessary because the accumulate function assumes x,y are relative to the center of the image, which for this process is the full image, not the stamp.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In fact, I think this is probably more complicated than necessary. You have a few things here that enforce this, which we could probably get rid of rather than do the following two lines.
The accumulate function has an option orig_center, which defaults to (0,0), but where you use imview.center. I think it's probably cleaner to not use that parameter at all and not shift the photon x,y by the center position. Then the accumulate function would use the normal x,y positions on the full image's normal grid.
So I recommend the following changes:

  1. Remove the next two lines and the above comment about them.
  2. Remove the orig_center option to accumulate below (given as an unnamed arg, now a kwarg -- it's your third arg).
  3. Remove the imview._shift call above. Which ironically means that the orig_center argument isn't doing anything, since the center of the view is 0,0, which is the default.
  4. Maybe even stop using the imview object entirely. I think that was just used in the old code because we wanted to shift the stamp centers to 0,0.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks, this sounds good. I'll make these changes and check I'm not breaking anything else.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

All looks good, and that's definitely made things simpler. Thanks for spotting it! Making the change now.

# Regular Sensors don't simulate this so don't accept recalc.
if imview.dtype in (np.float32, np.float64):
if isinstance(sensor, SiliconSensor):
sensor.accumulate(photons, imview, imview.center, resume=resume, recalc=recalc)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I already suggested removing imview.center here. But if I hadn't, then this should really be a named kwarg, not an arg, since otherwise it's not clear from the call what exactly it means.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks -- I'll go with removing imview.center as above.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And done with its removal.

RegisterImageType('LSST_PhotonPoolingImage', LSST_PhotonPoolingImageBuilder())


class PhotonPoolingTruthBuilder(TruthBuilder):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is this for? It doesn't seem to be used anywhere (tests or example config files).

It seems to be basically just recapitulating the regular TruthBuilder with a special case for incident_flux. I suspect we could do this another way that doesn't require so much duplicative code.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This (or rather something with the same effect) was a request of Jim's from above and discussed in the calls in January/February. With photon pooling there's no way of getting the realized_flux for photon objects as they go into the pool before being accumulated. This was an attempt to at least get the flux in all the photons for each object into the truth file, even if it includes some that isn't actually accumulated. Better than having 0. It's currently used in imsim-config-photon-pooling.py, which also disables the regular truth.

I could maybe try doing the sum of photons across the objects during the pooling loop in the image builder. But wouldn't that then need processImage instead of processStamp?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've got a new implementation which mostly uses the inherited code from the regular GalSim truth, only doing a little extra to sum the total flux in the photon shooting objects. I've also added a small assertion to the photon pooling test in test_image.py which checks that the nominal flux and incident flux are within 10% of one another. That could maybe be changed e.g. to within Poisson noise if you'd prefer.

if icrf_to_field is None:
icrf_to_field = create_test_icrf_to_field(boresight, det_name=det_name)
img_wcs = create_test_img_wcs(boresight, rottelpos)
img_wcs = create_test_img_wcs(boresight, rottelpos)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

duplicate

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks - removing the extra line now.

altitude=89.9 * degrees,
sky_pos=galsim.CelestialCoord(0.543 * galsim.radians, -0.174 * galsim.radians),
icrf_to_field=None,
rottelpos=np.pi / 3 * galsim.radians,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe more intuitive as 60 * galsim.degrees.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks, doing this now throughout the photon op tests.

icrf_to_field=None,
sky_pos=galsim.CelestialCoord(0.543 * galsim.radians, -0.174 * galsim.radians),
image_pos=galsim.PositionD(809.6510740536025, 3432.6477953336625),
stamp_center=galsim.PositionD(809.6510740536025, 3432.6477953336625),
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This doesn't make sense as a stamp center. An image pos can be arbitrary with lots of digits, but a stamp center is always either integers or half-integers. I guess it probably doesn't matter much for this test what the value is exactly, so maybe just change it to (809.5, 3432.5) to assuage my cognitive dissonance.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sounds good -- thanks! That definitely makes sense. Fixing this now.

"""Check the config interface to RubinDiffractionOptics."""

image_pos = galsim.PositionD(3076.4462608524213, 1566.4896702703757)
stamp_center = galsim.PositionD(3076.4462608524213, 1566.4896702703757)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ditto here and a few other places I guess.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixing throughout.

# Create a temporary ImageD to work in.
im1 = galsim.image.ImageD(bounds=imview.bounds)
if isinstance(sensor, SiliconSensor):
sensor.accumulate(photons, im1, imview.center, resume=resume, recalc=True)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think recalc and resume won't work here, since it's always a new image. We probably don't ever use this branch in imsim, but just in case, maybe issue a warning if resume is True or recalc is False that these don't work when the image type is not float32 or float64.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the suggestion -- adding warnings now.

@welucas2 welucas2 merged commit 8d9c804 into main Jun 5, 2025
4 checks passed
@welucas2
Copy link
Collaborator

welucas2 commented Jun 5, 2025

And done! Thanks everyone.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants