Help with CModel: returning all nans

jrmullaney · July 20, 2017, 3:09pm

Hi all,

We’ve been banging our heads against the wall with this one for a while, so could really use some help…

When we try to implement CModel for our obs_package (by mimicking the processCcd.py and cmodel.py files in obs_subaru) singleFrameDriver.py fails with the error:
RuntimeError: Unable to measure aperture correction for required algorithm 'modelfit_CModel_exp': only 0 sources, but require at least 2.

We’ve traced the cause of this error to CModel returning nans for all the sources it measures, including all the stars returned by the starSelector. As such, all stars are rejected when measureApCorr.py attempts to measure the aperture correction from the selected stars (around line 234).

Any suggestions as to why CModel may be doing this? We’ve checked the selected stars in our images, and they seem perfectly normal (i.e., not saturated, reasonable numbers of counts, nicely distributed around the image etc.).

If is helps, our obs_package is here: https://github.com/GOTO-OBS/obs_goto (apologies for the mess that is our processCcd.py config file). We’re using v13 of the stack, but we had exactly the same problem in v12.

Thanks in advance,

James

boutigny · July 20, 2017, 3:17pm

Have you set the following in your config file ?

config.charImage.measurement.plugins.names |= ["modelfit_DoubleShapeletPsfApprox", "modelfit_CModel"]```

jrmullaney · July 20, 2017, 3:34pm

Not exactly, but we have the following in our cmodel.py file:

import os
try:
    import lsst.meas.modelfit
    config.measurement.plugins.names |= ["modelfit_DoubleShapeletPsfApprox", "modelfit_CModel"]
    config.measurement.slots.modelFlux = 'modelfit_CModel'
    config.catalogCalculation.plugins['base_ClassificationExtendedness'].fluxRatio = 0.985

except (KeyError, ImportError):
    print("Cannot import lsst.meas.modelfit: disabling CModel measurements")

which gets called in processCcd.py with the line:

config.charImage.load(os.path.join(configDir, "cmodel.py"))

which is the same as that in obs_subaru.
My understanding is that this will do the same job as:

config.charImage.measurement.plugins.names |= ["modelfit_DoubleShapeletPsfApprox", "modelfit_CModel"]

in our processCcd.py config file.

price · July 20, 2017, 3:44pm

I don’t see anything obviously wrong in the config file.
Can you see which of the CModel flags are set in the catalog? That might give us a clue as to what the problem is.

jrmullaney · July 20, 2017, 4:11pm

Sure. This is what print(catalog.extract("modelfit_CModel_flag*")) returns (bullet points mine, obviously):

‘modelfit_CModel_flags_region_usedFootprintArea’: array([False, False, False, …, False, False, False], dtype=bool)
‘modelfit_CModel_flags_region_usedInitialEllipseMax’: array([False, False, False, …, False, False, False], dtype=bool),
‘modelfit_CModel_flag_noShape’: array([False, False, False, …, False, False, False], dtype=bool),
‘modelfit_CModel_flag_badCentroid’: array([False, False, False, …, False, False, False], dtype=bool)
‘modelfit_CModel_flag_apCorr’: array([False, False, False, …, False, False, False], dtype=bool)
‘modelfit_CModel_flags_region_usedInitialEllipseMin’: array([False, False, False, …, False, False, False], dtype=bool),
‘modelfit_CModel_flag_noShapeletPsf’: array([ True, True, True, …, True, True, True], dtype=bool)
‘modelfit_CModel_flags_region_usedPsfArea’: array([False, False, False, …, False, False, False]
dtype=bool),
‘modelfit_CModel_flag_region_maxArea’: array([False, False, False, …, False, False, False], dtype=bool),
‘modelfit_CModel_flags_smallShape’: array([False, False, False, …, False, False, False], dtype=bool),
‘modelfit_CModel_flag_region_maxBadPixelFraction’: array([False, False, False, …, False, False, False], dtype=bool),
‘modelfit_CModel_flag’: array([ True, True, True, …, True, True, True], dtype=bool)}

I’m just running it again to get the modelfit_CModel_{dev, initial and exp}* flags.

Thanks,

jrmullaney · July 20, 2017, 4:15pm

…and here are the rest:
for dev:

‘modelfit_CModel_dev_flag’: array([ True, True, True, …, True, True, True], dtype=bool)
‘modelfit_CModel_dev_flag_maxIter’: array([False, False, False, …, False, False, False], dtype=bool)
‘modelfit_CModel_dev_flag_trSmall’: array([False, False, False, …, False, False, False], dtype=bool)
‘modelfit_CModel_dev_flag_numericError’: array([False, False, False, …, False, False, False], dtype=bool)
‘modelfit_CModel_dev_flag_apCorr’: array([False, False, False, …, False, False, False], dtype=bool)}

for initial:

‘modelfit_CModel_initial_flag_numericError’: array([False, False, False, …, False, False, False], dtype=bool)
‘modelfit_CModel_initial_flag’: array([ True, True, True, …, True, True, True], dtype=bool),
‘modelfit_CModel_initial_flag_apCorr’: array([False, False, False, …, False, False, False], dtype=bool),
‘modelfit_CModel_initial_flag_maxIter’: array([False, False, False, …, False, False, False], dtype=bool),
‘modelfit_CModel_initial_flag_trSmall’: array([False, False, False, …, False, False, False], dtype=bool)}

for exp:

‘modelfit_CModel_exp_flag’: array([ True, True, True, …, True, True, True], dtype=bool)
‘modelfit_CModel_exp_flag_maxIter’: array([False, False, False, …, False, False, False], dtype=bool)
‘modelfit_CModel_exp_flag_trSmall’: array([False, False, False, …, False, False, False], dtype=bool)
‘modelfit_CModel_exp_flag_numericError’: array([False, False, False, …, False, False, False], dtype=bool)
‘modelfit_CModel_exp_flag_apCorr’: array([False, False, False, …, False, False, False], dtype=bool)}

Thanks,

Update: I’ve just done a quick uniqueness check on each of those flag arrays and can confirm that the arrays contain the same values throughout.

price · July 20, 2017, 5:07pm

So, looks like noShapeletPsf is the clue. The doc for that flag says: “the multishapelet fit to the PSF model did not succeed”.

So what’s your PSF model look like? Is your PSF photometry any good?

jrmullaney · July 20, 2017, 8:42pm

Thanks. It looks like we may be getting to the bottom of this…

Indeed, my PSF model looks a bit odd. I’d certainly expect it to be better defined than this (since there are some nice, clear, isolated stars in my input image) :

It returns PSF photometry (i.e., it gives reasonable-looking numbers), but I’ve not checked how good this is.

Any suggestions to help me improve my PSF model? I’ve previously tried using PSFEx, but this didn’t solve the problem (but am happy to try again).

Thanks,

price · July 20, 2017, 8:51pm

That looks rather undersampled. I don’t think any of our PSF models (especially not the default PcaPsf) handle undersampled data, and some important stages in our pipeline (warping) implicitly assume that we’re not undersampled.

jrmullaney · July 20, 2017, 9:05pm

Yes, I’m afraid our pixels are rather large: 1.2x1.2arcsec, and our seeing is likely to be around 2arcsec.

price · July 20, 2017, 9:22pm

So coadds aren’t going to be very good, nor PcaPsf. Try PsfexPsf — it has some (limited and possibly broken) support for undersampled data (you may need to turn on the oversampling feature).
CModel will need some tuning, if it’s going to work at all. Try looking at the modelfit_DoubleShapeletPsfApprox parameters (e.g., bring the inner radius down).

jrmullaney · July 20, 2017, 9:29pm

OK. Thanks very much for the suggestions. We’ll have a play around and hopefully we’ll be able to get it to work. I’ve switched PsfexPsf back on and it does, indeed, seem to give a better model PSF.

The telescope is a bit of a light-bucket, so we don’t need it to parameterise shapes well; the important part is reliable photometry. Hope this doesn’t end up being a complete deal-breaker for using the stack to process its data.

Thanks again, and do let us know if you think of anything else in the meantime,

J

jrmullaney · July 21, 2017, 8:54am

Just a quick follow up on the last post – not least as it may be useful to others working with undersampled PSFs…

In our processCcd.py we now have:

config.charImage.measurePsf.reserveFraction = 0.2
config.charImage.measurePsf.starSelector["objectSize"].sourceFluxField = 'base_PsfFlux_flux'
try:
    import lsst.meas.extensions.psfex.psfexPsfDeterminer
    config.charImage.measurePsf.psfDeterminer["psfex"].spatialOrder = 2
    config.charImage.measurePsf.psfDeterminer["psfex"].psfexBasis = 'PIXEL_AUTO'
    config.charImage.measurePsf.psfDeterminer["psfex"].samplingSize = 2.
    config.charImage.measurePsf.psfDeterminer["psfex"].kernelSize = 21
    config.charImage.measurePsf.psfDeterminer.name = "psfex"
except ImportError as e:
    print("WARNING: Unable to use psfex: %s" % e)
    config.charImage.measurePsf.psfDeterminer.name = "pca"

and in our cmodel.py we have:

import os
try:
    import lsst.meas.modelfit
    config.measurement.plugins.names |= ["modelfit_DoubleShapeletPsfApprox", "modelfit_CModel"]
    config.measurement.slots.modelFlux = 'modelfit_CModel'
    config.catalogCalculation.plugins['base_ClassificationExtendedness'].fluxRatio = 0.985
except (KeyError, ImportError):
        print("Cannot import lsst.meas.modelfit: disabling CModel measurements")

config.measurement.plugins["modelfit_DoubleShapeletPsfApprox"].minRadius = 0.1

This gives us a much better PSF model than shown in my previous post, but we’re still working on improving it. The good news is that CModel is now not producing all nans. It still does produce some nans, but there are sufficient finite values for the stack to “measure aperture correction for required algorithm”, thus solving the original problem.

Going forward, we’re going to tweak our psf fitting parameters to see if we can improve our model psf. We’ll also take a look at what can be adjusted in Cmodel.

Thanks for your help, @price - very much appreciated. You’ve saved us days of further headache!

price · July 21, 2017, 1:47pm

I’m very glad to hear that it’s moving forward again. Please let us know how it goes, and if you run into any further problems.

And while I’m thinking of it, you may want to grep the configs (e.g., the output of --show config) for lanczos (it’ll be combined with a number: lanczos3, lanczos5 or whatever) and override that with bilinear wherever you find it, especially:

<whatever>.warpAndPsfMatch.warp.warpingKernelName='bilinear'

That might salvage coadds for you. And I just noticed that we don’t have a config parameter controlling the choice of interpolation kernel for CoaddPsf so I just filed a ticket.

I’m not sure if changing the interpolation kernel is the right thing to do for the aperture measurements (which use this algorithm; perhaps @RHL can advise); for those you may want to set maxSincRadius = 0 so there’s no sinc interpolation.

jrmullaney · July 21, 2017, 3:43pm

Thank you. We’ll try your bilinear suggestion. And thanks for raising the ticket; I really appreciate you raising it, even though undersampled data won’t be a problem for LSST. Excellent community support!

For now, we’re running into some odd behaviour with the output of CModel. If we plot the positions of all the sources that CModel measures flux for (rather than returning a nan) on a warped image, we get the green points shown here(!):

(I’ve turned the contrast right down so it’s easy to see the positions of the measured fluxes, but rest assured, there are sources all over the image).

Clearly, there’s some weird positional dependence on whether CModel works or not. It’s going to be interesting getting to the bottom of this, but I feel we’re in a much better position than 24 hours ago!

I’ll keep you updated.

Thanks,

J

price · July 21, 2017, 3:49pm

I’m going to guess that that’s the PSF model again (its interpolation over the image). If so, it would be correlated with PsfMag - ApMag variations for stars.

I believe LSST has to worry about undersampling in the best seeing, but even if it doesn’t it’s something we certainly do worry about with HSC so this is a useful exercise for me too.

jrmullaney · July 24, 2017, 9:23am

Over the weekend I investigated the problem of the positional dependence of the PSF model. Unfortunately, I’m still unable to get modelfit_CModel_flux over an entire image, but I’ll summarise what I’ve found:

With characterizeImage.py (and using the config in the previous post) I am able to get modelfit_Cmodel_initial_flux for sources evenly distributed over the entire image, which is great. A small fraction (<5%) return nans, but these seem to be saturated sources, so I think that’s fine.

However, modelfit_Cmodel_flux returns large numbers of nans with a strong positional dependence (see below; red points are nans, green are finite fluxes):

On investigation, all the sources that have nan fluxes have modelfit_CModel_flag_region_maxBadPixelFraction flagged as True, whereas this is False for all sources with finite fluxes. It seems, therefore, that they are failing due some kind of bad pixel mask. However, nothing changes when I either set:

config.measurement.plugins['modelfit_CModel'].region.badMaskPlanes=[]
config.measurement.undeblended['modelfit_CModel'].region.badMaskPlanes=[]

or set

config.measurement.undeblended['modelfit_CModel'].region.maxBadPixelFraction
config.measurement.plugins['modelfit_CModel'].region.maxBadPixelFraction

to either 1, 2, 1000, or None.

I delved into meas_modelfit/include/lsst/meas/modelfit/CModel.h (There Be Dragons!) and came across this comment:

Following the initial fit, we also revisit the question of which pixels should be included in the fit (see CModelRegionControl)

which got me wondering whether modelfit_cmodel_initial is, for some reason, flagging whole regions of my images as bad. However, there doesn’t seem to be a (obvious) flag in the catalog schema that would indicate why this is happening.

Finally, adjusting the config.charImage.measurePsf.psfDeterminer["psfex"] parameters can help a little, but I always have large swathes of sources in my images that return nans for modelfit_cmodel_flux (and the relationship between parameter values and results seems rather chaotic).

So, my questions are:

Does my suspicion that modelfit_cmodel_initial is flagging many (contiguous) sources as bad prior to running model_cmodel_dev make sense?
Does anyone have any thoughts on why this would be the case?
Similarly, are there any flags that would tell me why modelfit_cmodel_initial is flagging them as bad?
Since we suspect this is related to the PSF model, is there some way I can extract or control how the PSF model is interpolated across the image?
Does anyone have any inspired suggestions as to how this could all be solved?
(I realise that last one is a wildly ambitious question!)

Thanks for helping me get this far, and thanks in advance for any suggestions.

jbosch · July 24, 2017, 3:30pm

The fit for modelfit_CModel_initial doesn’t flag pixels for later steps, but it does reset the pixel region we use in the fit and hence the places we might find bad pixels. Since you’ve set regions.badMaskPlanes to empty, the only way you can get bad pixels in the fit region is by having them fall off the edge of the image: themaxBadPixelFraction flag seems to be saying that the initial fits for these objects are so large that a large fraction of the pixels it wants to use are off the edge of the image. You could look at some of the modelfit_cmodel_initial_ellipse values to confirm that.

However, if your pixels are undersampled and/or your PSF model is poor, I think it may be pretty hard to get CModel working at all. In that regime, I think the only way you’ll get reliable photometry from the stack is with aperture photometry (and fairly large apertures).

price · July 25, 2017, 3:33am

The way I’d look at this would be to set a breakpoint in gdb (after compiling with scons opt=g) in this method and then poke around to see why maxBadPixelFraction is getting set to true (if you’re not comfortable with gdb, you could put a bunch of cerr << whatever lines). Check that the ellipse is sane, that footprint->getArea() == 0 is false, that _badPixelMask is 0 as you expect, that originalArea - spans->getArea() is 0 (because nothing got masked away).

jrmullaney · July 25, 2017, 11:11am

Thanks, both, for your suggestions.

Following up on the suggestion of @jbosch :
I’ve taken a look at the modelfit_CModel_initial_ellipse_* and modelfit_CModel_region_initial_ellipse_* values for all the sources in dmeRes.sourceCat (generated here). In most cases, there doesn’t seem to be any major differences in the modelfit_CModel_initial_ellipse_* values for the sources that have finite vs. nan values for modelfit_CModel_fit. For example, this is the distribution for modelfit_Cmodel_initial_ellipse_xx, where the red histogram is for nans, and the blue is for finite values (note, the histograms are normalised, in reality there are many more nans than finite values)):

The distributions are clearly different, but they span the same regions of parameter space. Does this mean that the initial fits are perhaps not too large?

There is a clearer difference in the distributions for the modelfit_CModel_region_initial_ellipse_* values, with sources with finite fluxes typically having 10x larger values for modelfit_CModel_region_initial_ellipse_{xx,yy} than sources with nans. There is, however, some overlap; see histogram of modelfit_CModel_region_initial_ellipse_xx, below:

I realise that this may not give much more insight, but I just wanted to give you an update on what I found. Maybe it’ll be useful for future reference if someone ever does want to tackle the undersampling problem.

I’ll take on board your comments about the current stack perhaps not being optimised for undersampled PSFs, although I’m too stubborn to give up this easily! As such, I’ll give @price 's suggestion a go, although delving into the C++ code is somewhat beyond my comfort zone.

Thanks again, your suggestions and patience are very much appreciated.