#7 HSC/LSST interface code

[msimet]

This pull request contains code that interfaces Stile with the HSC/LSST pipeline.

Most of the action in this PR happens in a new module for Stile called "lsst". It is found in stile/lsst. (This is actually kind of misnamed at the moment--it relies on HSC functionality that isn't in the LSST side yet.) There are two files of interest in that directory:

• sys_test_adapters.py, which makes objects that wrap the Stile tests in a way that makes them useful for LSST;
• base_tasks.py, which contains LSST-pipeline-like Task objects to pull data from the butler, coerce it into a format that Stile likes, and then call the wrapper objects from sys_test_adapters.py to run the tests. We've inherited from the LSST Task classes, which let us use the pipeline parser functionality to describe which data we're interested in: this is the 'visit ID' stuff from the pipeline documentation.

You actually run this code via some scripts in the (new) bin/ directory. The available scripts are:

• StileCCD.py : run CCD-level tests. This version assumes you're running on a CCD that has been used in a coadd and can use the coadd-level photometric calibration, etc. If there's no coadd associated with the CCD--if you only ran the single-frame processing--then you run:
• StileCCDNoTract.py
• StileField.py, StileFieldNoTract.py: run field-level (single-exposure) tests.

Four science choices that I wanted to highlight as needing particular input or ideas or verification from others:

• We have a "bright star selection" code block to pick bright stars. Right now it selects the top 10% of bright stars, but--due to some quirks in the way I coded this up--that will always run on a per-CCD level rather than on a per-field level even if we ask for fields. So it's probably better to have a hard cutoff. What should that cutoff be?
• One possible test is to run a galaxy-galaxy shear measurement around galaxies--not to do real science, and not weighted by redshift, just average shears. For that, we need a galaxy lens sample. Is there an obvious subset of galaxies we should be using for this?
• At the moment, when we need galaxy weights for correlation function-type tests, we just return 1 for all objects. That's probably a bad choice. Is there something in the catalog we should be using, or should we be pulling weights based on the shape measurement error, or something else?
• The flags that we use to exclude data (the flags in the list in 'removeFlaggedObjects') were chosen to be conservative right now--we're only looking at things we're fairly sure are well-measured. We will probably want to relax this in future (and make it configurable without editing the source code too).

This PR also includes some more systematics tests (shear around bright stars, star-star shape correlation function, galaxy-galaxy lensing type correlation function); they're all extensions of the existing CorrelationFunctionSysTest framework, so not much code there. I also added a rather unwieldy but flexible-to-the-user plotting function for the correlation function tests.

Still to be done:

• As Hironao noted in issue #22, we don't let the user select tests right now--the list is hard-coded.
• Tract and patch (coadd-processing) level tests. Shouldn't be too hard to code up, but Hironao has new code that is waiting on this branch.
• We don't generate randoms. We will need to eventually, but in the meantime we can't run any tests that use two-point correlation functions, because corr2 requires randoms for those.
• Two of the systematics tests included in sys_test_adapters.py are commented out of the default tests-to-do lists in base_tasks.py. That's because of issue #23, so they should be included once that's finished.
• At the moment we don't have any unit tests. I'm not sure exactly how we do that without assuming some available data (or, really, even once we have some data to run it on)--perhaps @TallJimbo or @HironaoMiyatake could advise and I could add that to this branch?

I want to especially thank @TallJimbo and @HironaoMiyatake who put a lot of work and time into this PR as well, which may not be represented in the git history for these files.

[TallJimbo]

I do plan to look at this at some point but likely won't have a chance until later this week (Thursday or Friday). I'm traveling and my laptop power cable has died, so I'm quite literally working on borrowed time until a replacement arrives.

[msimet]

Yikes, Jim! Hope it lasts as long as you need it. In any case, I'm happy to take comments whenever you have them.

[rmandelb]

It is found in stile/lsst.

I recommend renaming to stile/hsc since it's a misleading name. If that functionality goes to the LSST side, we can always rename again.

So it's probably better to have a hard cutoff. What should that cutoff be?

Are you suggesting a hard cutoff in S/N, in maximum number of stars, or ...?

Is there an obvious subset of galaxies we should be using for this?

That's tricky. Most of the obvious samples I can think of are low density enough that at the CCD level, you might have an average of <1 per CCD.

At the moment, when we need galaxy weights for correlation function-type tests, we just return 1 for all objects.

We should weight based on the shape measurement error and shape noise added in quadrature. This is not too far from uniform weighting so it's not like what you're doing right now is horribly wrong. But probably better to update this.

The flags that we use to exclude data (the flags in the list in 'removeFlaggedObjects') were chosen to be conservative right now--we're only looking at things we're fairly sure are well-measured. We will probably want to relax this in future (and make it configurable without editing the source code too).

Open another issue for this? And for the rest of the still-to-be-done list?

As Hironao noted in issue #22, we don't let the user select tests right now--the list is hard-coded.

This is no longer true, right?

in the meantime we can't run any tests that use two-point correlation functions, because corr2 requires randoms for those.

Do you mean two-point correlations that involve a galaxy density? You are doing shape two-point correlation functions, I thought.

Unit tests: would be good to get some very small amount of default data that could be used for this.

[msimet]

I recommend renaming to stile/hsc since it's a misleading name. If that functionality goes to the LSST side, we can always rename again.

Done.

So it's probably better to have a hard cutoff. What should that cutoff be?

Are you suggesting a hard cutoff in S/N, in maximum number of stars, or ...?

I'm open to suggestion, basically. (I think S/N would make more sense, probably, or raw flux, rather than number of objects.)

That's tricky. Most of the obvious samples I can think of are low density enough that at the CCD level, you might have an average of <1 per CCD.

We could always just try a random selection, or semi-random above a flux cutoff--I think the important thing is to have some real objects for the test, but using all the galaxies could take too long.

We should weight based on the shape measurement error and shape noise added in quadrature. This is not too far from uniform weighting so it's not like what you're doing right now is horribly wrong. But probably better to update this.

Right, will do later today probably. (These shapes are distortions, yes?)

Open another issue for this? And for the rest of the still-to-be-done list?

Okay, later today as well.

As Hironao noted in issue #22, we don't let the user select tests right now--the list is hard-coded.

This is no longer true, right?

Right, yes: for those not watching the other issue, you can change the tests now--there are instructions if you pass --help to the scripts.

in the meantime we can't run any tests that use two-point correlation functions, because corr2 requires randoms for those.

Do you mean two-point correlations that involve a galaxy density? You are doing shape two-point correlation functions, I thought.

Sorry, yes, two-point as in position-position. I think position-shear will work without randoms but can't be sure till I fix issue #23.

Unit tests: would be good to get some very small amount of default data that could be used for this.

I'm curious what the LSST folks are doing already to test their code...

[HironaoMiyatake]

We should weight based on the shape measurement error and shape noise added in quadrature. This is not too far from uniform weighting so it's not like what you're doing right now is horribly wrong. But probably better to update this.

Right, will do later today probably. (These shapes are distortions, yes?)

My crude weight used for my tests so far is wt = 1./(0.362 + eerr2), where I think 0.36 is derived from SDSS in some of Reina's paper (@rmandelb could clarify). eerr is error in distortion, 1 component.

Unit tests: would be good to get some very small amount of default data that could be used for this.

The problem is that the HSC data is not public. Now is it time to set up a private repository? Or we could use somewhere else which is protected by a password.

[rmandelb]

So it's probably better to have a hard cutoff. What should that cutoff be? Are you suggesting a hard cutoff in S/N, in maximum number of stars, or ...? I'm open to suggestion, basically. (I think S/N would make more sense, probably, or raw flux, rather than number of objects.)

I don't think raw flux makes sense. S/N is probably best.

That's tricky. Most of the obvious samples I can think of are low density enough that at the CCD level, you might have an average of <1 per CCD. We could always just try a random selection, or semi-random above a flux cutoff--I think the important thing is to have some real objects for the test, but using all the galaxies could take too long.

Okay, perhaps a random subsample for those above a S/N cutoff.

We should weight based on the shape measurement error and shape noise added in quadrature. This is not too far from uniform weighting so it's not like what you're doing right now is horribly wrong. But probably better to update this. Right, will do later today probably. (These shapes are distortions, yes?)

Yes, distortions.

[rmandelb]

Regarding a small amount of data: does it have to be HSC data? If the formats used for outputs is the same for some public S-Cam data, we could run the pipeline on some S-Cam data (like one CCD) and put that output in the repo for testing.

[HironaoMiyatake]

I think SC data is okay. The ACT cluster data is public, so we can use it. However, I have to run soon for a short trip during this Friday and weekend, so it might be better to make a new issue for it, depending on time scale of this PR.

[rmandelb]

I think it’s okay to open a new issue and do it next week or whatever. It’s fine as long as we have a plan!

[HironaoMiyatake]

Okay, I'll open an issue!

[TallJimbo]

Review complete; all my comments are in the diff comments.

[msimet]

My crude weight used for my tests so far is wt = 1./(0.362 + eerr2), where I think 0.36 is derived from SDSS in some of Reina's paper (@rmandelb could clarify). eerr is error in distortion, 1 component.

Yes, 0.36 was what they measured in Reina's paper (0.365 actually I believe). However, that was for the sample after some cuts--I've put it in for now (with max(e1_err,e2_err) as eerr) but we might have to think more carefully about it in the future.

[msimet]

Noting for those not looking at the diff that what I called StileField.py and StileFieldNoTract.py are now StileVisit.py and StileVisitNoTract.py thanks to some points @TallJimbo made above about terminology.

[msimet]

Does anyone else have any comments on this PR? Does anyone want to look over the code that's been added due to comments? If nobody speaks up, I will merge tomorrow evening.

[msimet]

Pursuant to some discussion on branch #18, I just changed the call signature for the SysTestAdapters here to include the config objects from the Task.

[msimet]

I just realized we needed slight fixes to the tests and example scripts here, so that's done. New plan: merge tomorrow.