mdschneider / jif Goto Github PK

When running an MCMC chain with the Spergel profile in a bulge+disk model, I frequently get RuntimeErrors from the galsim.drawImage() method because the requested FFT grid size is too large. This is because certain values of nu require very large grids to render.

I don't know what the right solution is here. At the moment, the code catches the RuntimeError and returns a very small log-likelihood as if these values of nu are just invalid parameter values. But, is this really how we want to treat these models?

When creating space and ground postage stamps of local groups do they need to have:

1. the same pixel scale?
2. the exact same RA, Dec extents?
3. the same WCS-Images coordinate orientation (or is an associated WCS header sufficient)?

Currently the PSF is output as fits and png files but we should also include this information in the hdf5 data structure.

This may be a non-issue but currently galsim_galaxy appears to generated a single noise value rather than a pixelated noise map (as is done for the real data), or at least just a single value is what is stored in the hdf5 file. Is Roaster is capable of distinguishing between the two types of input? Or do I need to add this capability (for now it could be as easy as reducing the real image noise map to a single value)?

@mdschneider,
I got the group identification and postage-stamp creation scheme working on real HST data. We just need to discuss how to formalize the output data structure. I was thinking it should be something like:

• group (perhaps have two branches after this for Space and Ground)
  • group (x,y) pixel bounds
  • perhaps pixel data postage stamp for group
  • group membership (i.e. objects that should be simultaneously analyzed)
    • respective member properties (e.g. coordinate, size, shape, etc.)

Let me know if you have any input on the matter.

Both the noise model and the zeropoint require exposure times. These are currently set independently and could easily be made inconsistent.

We don't want to entirely reinvent the wheel with our data structure. I think the only astronomical survey currently using the HDF5 format is LOFAR

Plus there may be similar/existing data structures and tools that we want to be compatible with. From what I know and given a brief search there are a number of examples:

1. Perhaps we should have a dataset format that is easily comparable with python pandas. When you create an hdf5 file from within pandas it conveniently creates column names either based on user input or column names already associated with the pandas data frame. This does not appear to be so convenient for h5py or hdf5 in general.

• Create a simple hdf5 file in pandas with column names and look how the hdf5 file is structured to see if this structure can easily be incorporated in what we are doing.

2. Similarly PyTables is something to investigate.

3. NASA Goddard has been using hdf5 datasets for representing some of their geo-gridded data related to their earth sciences programs. They even have a GUI interface call Panoply which has potential to be very useful for viewing our astronomical data.

4. Another NASA project using and visualizing HDF5 files is the Space Weather Explorer.

It is good that NASA is already using this on a few projects since we will want WFIRST to use this format.

5. One of the better graphical interface tools is Paraview.

There is a list of many of the HDF related software here. Note that I know this list is incomplete.

Since the whole field fits files are saved in either the Observation or Sextractor group it would be nice if we could just use region references in the Group group to define the group postage stamps (clearly need to come up with a naming convention that does not use "group"). This could potentially cut the storage size in ~half.

One key disadvantage in doing so is that the returned data will be a flattened 1D array.

sheller_great3.py uses a kludgey noise model right now for setting up GREAT3 galaxies to be processed by Roaster. We need to do something more accurate.

Roaster has a parameter optimization step prior to MCMC to attempt to reduce the burn-in period. But this optimization step often fails. We need a reliable likelihood optimizer in Roaster that works for a range of input galaxy images with both isolated and blended sources.

Currently there are attributes that specify filter properties (e.g. central, FWHM; see c86ece1) but we should really just have another subgroup that specifics filter throughput curves and just have an observation specific attribute referencing a specific filter name.

Need to add a function or two to RoasterInspector that will generate model vs data plots for a given segment. These can be a simple four (8 if space & ground) panel plot with: data, model, model+noise, residual.

Some of this code exists in showseg.py.

The galsim_galaxy.py module contains our galaxy image model. If run as a command-line script it generates some test data for us to look at and use as simple inputs to Roaster.py.

The output HDF5 file saved in galsim_galaxy.py here needs to be updated to match the actual formats that @wadawson is designing for handling blended groups with multiple epochs/instruments.

We can then continue to use this functionality to make Roaster inputs that are generated with the same model as used to fit the data - the simplest possible test in some ways.

With the multiple importance sampling approach for combining epochs and passbands we have less urgent need for chromatic modeling in GalsimGalaxy. Separating the simpler achromatic implementation should allow for easier testing as well.

This is needed as input for for some of the DM Stack shape measurement routines. For now we will assume perfect knowledge of the PSF. Once we move to analyzing the GREAT3 data we can account for imperfect knowledge of the PSF.

Compound types are more efficient than generic multidimensional arrays when you have repeated structures in you dimensions (e.g. N parameters with M MCMC samples each). Only the "columns" of the dataset that are referenced in any particular call are reading (rather than the whole dataset) so it is faster as well. Not to mention compound arrays are also more friendly to use.

See for example how sheller.py stores the stamp_objprops datasets.

The LSST effective PSF (i.e. PSF as realized on the image pixel grid) currently output by galsim_galaxy.py (commit: fdc618a) seems a bit under sampled, see the following log-scale image:

Without the log-scale only the four central pixels light up. This is probably an indication that the combination of LSST optics/atm/detector inputs should be tweaked.

We want to model blended images as seen in different wavebands. We can use the galsim.Chromatic functionality for this, but this needs to be implemented in galsim_galaxy.py.

The more difficult issue is how to parameterize the galaxy SED so we can sample SED model parameters the same way we do the morphology and orientation. One simple idea is to make an SED model that's a linear combination of SED templates (or maybe SED PC modes etc.). These parameters should also be a part of galsim_galaxy.py.

While investigating checker-cross pattern artifacts in the images I found that galsim_galaxy.py is currently generating inconsistent LSST and WFIRST images. As you can see from the image below which shows a screen capture from ds9 where I have used ds9 Frame>match>WCS to match the LSST image to the WFIRST image. When doing this they should be the same size. Also note that the WFIRST image is 416x416 pixels where as the LSST is 28x28 pixels.

This might be why we are getting counterintuitive posteriors when running the roaster on the two images.

To test our ideas for discretizing the centroid values implemented in, e.g., Spergelets, we need MCMC sampling that can handle discrete parameter values. emcee cannot do this. We could use PyMC perhaps. But maybe CosmoSIS would be even better given the existing wrappers to many samplers?

Have it write out a file in the same format as that in the preroaster branch.

The objective is to plot the postage stamps for individual segments from space and ground side by side. Really just for ease of visual inspection.

There is an error (about pickling) when the --nthreads argument to Roaster.py is greater than 1.

@mdschneider , I just wanted to touch base and see how complex you want the current hdf5 data to be for input into the roaster. If you have been following along the preroaster branch you can see that I am currently convolving the HST image with a "ground" Gaussian seeing kernel to define the grouped objects. This avoids the WCS/pixel/orientation complexity in defining postage stamps in the various images. I was thinking that such simplifications would be nice for you working on this data with the roaster but I thought I would check with you first.

What is the bare minimum you would like to have for our first analysis of the data with the roaster? I am imaging:

• "ground" postage stamps (image-background, noise, group pixel mask)
• space postage stamps (image-background, noise, group pixel mask)
• corresponding coordinates of the objects detected in the space image

Note that I can give you WCS header info but this is not necessary since the "ground" and space stamps have the same pixelation/orientation/WCS. I have also run sextractor on both "ground" and space images so I can give you what ever properties you would like for each object.

A first step for our DESC PSF Task Force is to make a simple numerical example of the marginal ellipticity distribution of a galaxy image when a parametric PSF is marginalized.

So far, we have only allowed for fixed (or known) PSFs. We can start by marginalizing just the isotropic PSF FWHM and ellipticity, and then evaluate whether we need more complicated PSF models for the next step. But this should get the code functionality in place.

The galaxy parameters and selection of galaxy types is incomplete.

This information is needed as part of the DM stack analysis to create the variance image of the MaskedImage data structure. We could take a variance image as input but for now I think it is easiest just to create the variance image as the Image/gain.