Yusra, do you have time for a review of some the shapelet PSF-fitting code? I finally tracked down the bug in the shapelet PSF fitting code you were working on a few months ago - the model evaluation was just wrong, as I wasn't incrementing a counter in the amplitude direction when evaluating the different shapelet terms. Aside from that bug, the only change I've made in the code you were working with is to merge the two unit tests for PsfFitter class, modify the test you wrote a bit, and adapt to changes on DM-1188, which affected the low-level model evaluation code in the shapelet package (that's what helped me find the bug).

There's a GitHub Pull Request linked from this page that will hopefully make the review easier.

If you would have a hard time getting to this in the next week or so, just let me know. I can find another review, but I wanted to give you the first shot at it since you did spend some time on this code already.

That's exciting that you found the bug. Before I read through the code, I found something while running the unit test:
data/psfs/great3-1-fits doesn't match the pattern *.fits so this image wasn't being tested. The test fails on this one, but passes if you increase the tolerances["ellipse"] from 1E-3 to 1E-2. Whereas we don't want all 200 great3 psfs in the unit tests, it'd be good to check that PsfFitter.apply() can run on all of them. I see PsfFitter.evaluate() creating some images with NaNs. I put all 200 psf thumbnails of the same size and naming structure in '/lsst8/yusra/test_shapelet/psfs') on lsst-dev so if you change:
DATA_DIR = os.path.join('/lsst8/yusra/test_shapelet/') in testPsfFitter.py the unit test will run on all 200 images.

Unrelated to this ticket:

in testTruncatedGaussian.py I get a

 return scipy.integrate.dblquad(func, 0.0, scipy.integrate.Inf,

AttributeError: 'module' object has no attribute 'Inf'

I really thought these were already changed to numpy.Inf in a previous ticket...

A unit test in shapelet had a tolerance that was too low:

[yusra@lsst-dev .tests]$ more testMatrixBuilder.py.failed 

tests/testMatrixBuilder.py

...F...

======================================================================

FAIL: testRemappedConvolvedShapeletMatrixBuilder (__main__.MatrixBuilderTestCase)

----------------------------------------------------------------------

Traceback (most recent call last):

  File "tests/testMatrixBuilder.py", line 220, in testRemappedConvolvedShapeletMatrixBuilder

    self.assertClose(numpy.dot(matrixD, coefficients), checkVector, rtol=1E-13)

  File "/lsst/home/lsstsw/stack/Linux64/utils/9.2+2/python/lsst/utils/tests.py", line 328, in assertClose

    testCase.assertFalse(failed, msg="\n".join(msg))

AssertionError: 1/50 elements differ with rtol=1e-13, atol=2.22044604925e-16

0.0156680409512 != 0.0156680409512 (diff=1.57165946923e-15/0.0156680409512=1.00309890313e-13)

I'll post comments on the code soon.

I just pushed fixes for the scipy.interpolate.Inf issue and the test data filename issue.

I couldn't reproduce the tolerance problem in shapelet, but I'll fix it anyway (likely on DM-1188), as I'm probably just getting different random numbers than you for a reason that doesn't matter (e.g. different NumPy versions or platforms).

I'll take a look at the other PSF images you put together. I certainly agree that this should be running with no NaNs and no exceptions on all of these, even if the fits aren't necessarily very good. Thanks!

I opted to write line level comments here rather than in github:

include/lsst/meas/multifit/psf.h:

• The docstrings for the PsfFitterComponentControl are split: See http://lsst-web.ncsa.illinois.edu/~yusra/meas_multifit/classlsst_1_1meas_1_1multifit_1_1_psf_fitter_component_control.html

src/psf.cc:

• I couldn't find where PsfFitterPrior is used by the code anywhere. Is to possible to add a test that uses it?

• I'm assuming PsfFitterModel and PsfFitterPrior are defined in psf.cc because PsfFitter is the only class that will use ever them.

tests/testPsfFitter.py:

• line 256: The comment is an old note from me to you. It could go.
• line 260: modelImageF is not used. (Note for future: it'd be nice if PsfFitter.apply could run on images with double pixels. – The reason this test needed a double and a float version of the psf image is because lsst::shapelet::MultiShapeletFunctionEvaluator::addToImage() takes doubles and psfFitter.apply() takes lsst::multishapelet::Pixel = floats. )

The docstrings for the PsfFitterComponentControl are split: See http://lsst-web.ncsa.illinois.edu/~yusra/meas_multifit/classlsst_1_1meas_1_1multifit_1_1_psf_fitter_component_control.html

I don't see what's wrong here. If you just mean that all the members appear twice on the page, once at the top, and then again at the bottom, I think that's an unavoidable thing Doxygen does with all classes. It's just more obvious here because the "brief" documentation is all any of these have.

I couldn't find where PsfFitterPrior is used by the code anywhere. Is to possible to add a test that uses it?

Good catch! I'd intended to construct one and put it in PsfFitter::_prior when any of the priorSigma parameters on the control object were nonzero and finite, but it looks like I'd forgotten that fairly important step. I'll fix that, and make sure it all still works.

I'm assuming PsfFitterModel and PsfFitterPrior are defined in psf.cc because PsfFitter is the only class that will use ever them.

Yes. And the fact that their public interface is purely defined by their base class, so as long as there's a public way to construct them, the user doesn't ever need to see the class definition.

line 256: The comment is an old note from me to you. It could go.

Will do.

line 260: modelImageF is not used. (Note for future: it'd be nice if PsfFitter.apply could run on images with double pixels. – The reason this test needed a double and a float version of the psf image is because lsst::shapelet::MultiShapeletFunctionEvaluator::addToImage() takes doubles and psfFitter.apply() takes lsst::multishapelet::Pixel = floats.

Agreed. The real limitation here is that the base Likelihood class only takes single precision floats, so the underlying work all has to be done in single precision until that's addressed. But I think a convenience method to do the conversion would probably be valuable here.

I created a new example script to make running on lots of PSFs at once easier, and I've now run on all 200 GREAT3 control/ground/constant PSFs, and with a small (but important) fix to the optimizer all of them now run to completion.

It's very slow at high shapelet order, which should be addressed in the future with DM-1216. In the meantime, I've added a warning about high shapelet order in the docs for the control object and modified the defaults to fit an elliptical double-Gaussian, which (from the very superficial inspection of the GREAT3 runs) seems like a decent compromise.

I also rebased the branch on master, to stay compatible with the latest afw and Russell's xy0 changes there, so be careful if you pull.

Yusra AlSayyad, anything else to add, or is this ready to merge to master?

OK. Nothing else to add; ready to merge.

What I meant "about the docstrings for PsfFitterComponentControl is ", "that they show up as double quote", "delimited chunks like this." Look at radiusFactor for example.

Ah, right. That's a known issue with documentation for any attributes defined with LSST_CONTROL_FIELD (I think there's a trac ticket for it, but maybe not a JIRA one).

Thanks!