jeff-regier / celeste.jl Goto Github PK

A patch is the region of images affected by a source. The code currently allows the user to specify the size of the patch for each source. But there isn't code yet that picks a good size patch, based on estimates of the brightness and size of the source.

Since in real images the sky background and calibration vector vary across the image, perhaps these quantities should be matrix-valued, not scalar-valued, in ElboDeriv.

I think iota should definitely be matrix-valued, since it varies quite a lot across the sky. We can pass around a CoordInterpGrid object instead of the full matrix, since that's all that's saved in the FITS file anyway.

The calibration vector -- which is the only variable part of epsilon -- only varies across a row, and it doesn't vary much (see attached). Perhaps we could get away with using its average if you wanted to.

Jeff, any thoughts before I jump in?

Polygons.jl is for finding the intersection of a circle with an arbitrary polygon. That's how you need to match sources to tiles if you do it in world coordinates.

Note that the WCS mapping takes world circles to pixel coordinate ellipses, but that tiles are pixel coordinate tiled squares. It will be much easier to get a conservative mapping from source to tile by mapping a world circle to a conservative pixel circle by just making the radius of the pixel circle the major axis of the ellipse and then finding its intersection with the square tiles.

In addition to this, we can speed up the ELBO quite a lot by making NaN any pixels that receive little light based on a rendering of the catalog initialization. The tile to source mapping is only then necessary for choosing the extent of the initial rendering, so we don't have to pay the price of being conservative for the whole ELBO.

Document the inputs and outputs of the core functions in Celeste.

...and replace it with the new trust-region-based method.

Benefits:

• Fixes the unit tests, which are failing now on Travis-CI because of some random NLopt problem
• Reuses the existing unit tests to validate the new Hessian-based way of optimizing
• I can finish #16 pretty quickly once we can set D > 2. The model should work much better once color doesn't misfit the data so badly.
• I can finish #14 pretty quickly once we can set D > 2, by making the color model for stars rich enough to represent quasars too. (Effectively "stars" come to represent all point sources.)

In order to read real frames which sometimes have bad pixels, we need to be able to assign zero weights to certain pixels in ElboDeriv.

Through pre-processing, transform the SDSS dataset into "tasks", 1 per known astronomical object. The tasks are bundled into JLD files. Each bulkan/brick gets its own JLD file. For each astronomical object, the corresponding task contains everything necessary for processing that object, including

• the tiles/pixels relevant to it
• the parameters for this astronomical object from an existing SDSS catalog (for initialization)
• the catalog entries for all neighboring astronomical objects (treated as fixed, for now)
• the point spread functions for the tiles
• the gain and calibration constants for the tiles
• local linearized world coordinates for each tile / subimage

For large galaxies, we probably want downsampled pixels rather than the original pixels.

Neighboring astronomical objects are those that may contribute photons to any of the included tiles.

In the future, we may want to include just identifiers for the neighboring astronomical objects, rather than the SDSS catalog entries for those neighbors. That would allow for an iterative procedure, where better estimates of an object's neighbors' parameters also improve estimates of the object's parameters.

The PSF varies across an SDSS field, and should be represented somehow (e.g., each source could have its own PSF learned at its catalog location).

This will make it easier to create tiled blobs that contain only the data needed for a particular source and save passing around a lot of data.

It won't compile easily on NERSC, and it's less accurate than ForwardDiff.jl anyway.

Using exact Hessians, recreate the results from the ICML submission for the test dataset of 51x51-pixel stamps.

Looks like bin/benchmark.jl is taking 6x as long now to compute the elbo as it did at the end of May. Deepcopy operations may be to blame:

  9909 deepcopy.jl                _deepcopy_t                                36
  9999 ...ste/src/CelesteTypes.jl +                                         619
 10171 ...eleste/src/ElboDeriv.jl accum_pixel_source_stats!                 558
 10248 ...eleste/src/ElboDeriv.jl accum_pixel_source_stats!                 567
 10668 deepcopy.jl                deepcopy_internal                          24
 11029 ...eleste/src/ModelInit.jl initialize_model_params                   394
 11043 ...eleste/src/ModelInit.jl initialize_celeste                        347
 11099 ...leste/src/SampleData.jl gen_n_body_dataset                        159
 13327 ...eleste/src/ElboDeriv.jl get_tile_sf                               813
 14513 ...eleste/src/ElboDeriv.jl accum_pixel_source_stats!                 569
 14685 ...ste.jl/bin/benchmark.jl small_image_profile                        15
 15183 ./boot.jl                  include                                   245
 15183 ./client.jl                process_options                           285
 15183 ./client.jl                _start                                    354
 15183 loading.jl                 include_from_node1                        128
 15183 profile.jl                 anonymous                                  14
 15934 ...eleste/src/ElboDeriv.jl accum_pixel_source_stats!                 570
 23063 ...eleste/src/ElboDeriv.jl accum_galaxy_pos!                         453
 26882 ...eleste/src/ElboDeriv.jl accum_galaxy_pos!                         454
 69416 ...eleste/src/ElboDeriv.jl accum_pixel_source_stats!                 518
155984 ...eleste/src/ElboDeriv.jl expected_pixel_brightness!                650
156413 ...eleste/src/ElboDeriv.jl tile_likelihood!                          715
165874 ...eleste/src/ElboDeriv.jl get_tile_sf                               815
195042 no file                    anonymous                                 820
195045 multi.jl                   anonymous                                1279
195045 multi.jl                   run_work_thunk                            621
195046 multi.jl                   run_work_thunk                            630
195046 task.jl                    anonymous                                   6

According to spectrographs, the objects at

• 5.1439, 0.0885
• 5.5156, 0.1048
• 5.7299, 0.0265, and
• 5.5853, 0.1438

are red dwarfs in our collection of stamps that Celeste incorrectly labels galaxies. Presumably our current color prior---a 2 component mixture model---doesn't have enough components to adequately model stars outside of the main sequence (e.g. red dwarfs).

http://skyserver.sdss.org/dr12/en/get/SpecById.ashx?id=1258871587503892480
http://skyserver.sdss.org/dr12/en/get/SpecById.ashx?id=1261149226344146944
http://skyserver.sdss.org/dr12/en/get/SpecById.ashx?id=1261148676588333056
http://skyserver.sdss.org/dr12/en/get/SpecById.ashx?id=1258898800416679936

I came across an NLOpt failure using the free transform and spend a lot of the afternoon tracking it down with autodiff. I managed to find a discrepancy between the hand-coded and forward autodiff derivatives. The discrepancy was actually due to a numerical error in the autodiff, not in in the hand-coded derivatives. However, the source of this error could manifest in subtle ways in the numerical optimization routine and we should keep it in mind.

I isolated the problem to the gen_gamma_kl function. When r1 and r2 are very large, this function involves differences of very large floating point quantities that need to cancel each other out, which is numerically unstable. This turns out to affect forward autodiff more than the hand-coded derivatives, though of course it potentially makes the actual value unstable as well.

Normally this isn't too bad because the box constraints keep r1 and r2 from being very large. However, due to having re-parameterized the brightness in terms of its mean and variance, the box constraints on the transformed parameters were no longer constraining the actual gamma shape and scale parameters.

I've returned to simply using positivity constraints on the brightness parameters and the problem has gone away for now. But if you do manage to find a reproducible instance of an NLOpt failure, let me know, because now we have some more powerful tools for diagnosing them.

Some model quantities are currently hard-coded (e.g. the number of PSF and galaxy components). Others have global constants defined but hard-coded values are sometimes still used (e.g. the mu values in the VB parameters object). With two of us making changes, it'll be less error prone to go through the code and encode all these as global constants.

All the test images in Synthetic.jl use either an identity WCS or a symmetric one. For testing purposes, we should replace them all with an actual asymmetric WCS.

To run Celeste.jl on the whole SDSS dataset, we need a function that each task can call to fetch just the relevant parts of the dataset. As input, this function takes a task id and I guess the total number of tasks. It returns

1. a subset of the astronomical objects (stars and galaxies) in existing catalogs for this task to optimize,
2. all the "tiles'' (a.k.a. "super-pixels"---subimages that are 10 pixels wide and 10 pixels high) that are "near" any astronomical objects in the subset, and
3. any additional astronomical objects that are near any of these tiles. (The task isn't responsible for optimizing these astronomical objects, but it may use them to help explain the data, i.e. the tiles.)

Encode source locations in WCS coordinates, not pixel coordinates. Stamps are pixel aligned, so pixel coordinates were fine before. Frames in different bands aren't exactly pixel aligned.

The SDSS pipeline doesn't approximate the PSF exclusively with Gaussians, so we need to fit a Gaussian mixture model ourselves. GaussianMixtures.jl looks useful.

with reparameterizations, rather than projection

Not just variational parameters now

The test test_elbo_invariance_to_a() actually contains an extremely particular set of parameters that allow the test to pass, and the VB solution is not at all invariant to a in general. For example, changing the fluxes in the unit test to 1000 rather than 100 times the sample fluxes causes test_elbo_invariance_to_a() to fail.

For an even more dramatic failure, try fitting only the star parameters to a dataset generated by gen_sample_star_dataset(perturb=false) with a fixed at 0.01, 0.5, and 0.99. Here are the fit star brightnesses:

  a = 0.01  a = 0.5  a=0.99
  4470.74   5304.13  2422.03    
  1581.14    987.354    0.109955
  2558.08    393.688    0.591057
  2646.56      0.845658 0.79425 
 18480.3   18627.8  17495.4  

I've explored a few other examples and it seems like a dependence is more common than independence. I'm not sure exactly what's going on. I think you're right that the coupling must be in the lower bound term. I can imagine this causing problems for optimization.

Change the constrained parameterization to an unconstrained one and update the analytic derivatives.

Modeling the PSF in Celeste with rasterized "eigenimages" rather than a mixture-of-Gaussians. Use cubic or lanzos interpolation.

Advantages:

• More accurate---and we've heard from Dustin via David that the mixture of Gaussians doesn't fit well.
• Faster to convert from the existing SDSS PSF format to a rasterized image than to fit a mixture of Gaussians.

At first, use the eigenimages (RawPSF) from the SDSS pipeline. For DECaLS, we'll have to learn these rasterized images ourselves, probably from within the model rather than from eigen decompositions.

We've still have NLopt in the REQUIRE file, but presumably it would be easier to deploy our code on NERSC with fewer dependencies, and installing NLopt.jl takes compiling some C code.

@rgiordan , you were saying we still need NLopt for a couple of unit tests, right? Are those tests essential? If this is a problem with negative curvature, might the trust region method have the same problems on the actual images that it has in the unit tests?

As dynamic objects they contain too much information that needs to be internally consistent, and deviations are hard to debug.

There are at least a couple people in the AMPLab who want to use some of this functionality.

I'll fix this, I'm in there anyway.

Test Celeste on all of stripe 82, not just on stamps from it.

rather than a gamma distribution

We should decide how to compute exact Hessians. Hessians may be computed either 1) manually, 2) automatically, by JuMP/ReverseDiffSparse, or 3) automatically, by forking JuMP and modifying it to parse a BUGS-like specification of our model. To decide, it'd helpful to understand whether our objective function can be written in the current JuMP modeling language.

Jeff -- what do you think about just writing a simple python wrapper around tractor to read in the existing catalog and write to a Celeste-readable csv or something?

You can take a look at what it does in tractor/sdss.py/_get_sources. It mostly looks like bookkeeping. This would also let us deal cleanly later with alternative initializations.

The objects at the following ra-dec coordinates are quasars in our dataset of "stamps".

• 5.0820, 0.9728
• 5.0075, 0.9102
• 5.2982, 0.0274
• 5.6378, 0.1212
• 5.0161, 0.1945
• 5.1439, 0.0885

I added a field h to SensitiveFloat (defined in CelesteTypes.jl). It has the same dimension as the first derivative d. We can store the 2nd derivatives (i.e. the diagonal of the Hessian) there.

For testing, I suggest modifying test/test_derivs.jl to also test these new 2nd derivatives. We can either extend test_by_finite_differences to also test second derivatives, or make a new function, e.g. test_second_deriv_by_finite_differences. In either case, we'd want to use GSL.deriv_central again (as we've done in test_by_finite_differences), but this time on the derivative values (the SensitiveFloat.d values, that is, which we've already verified) rather than the objective.

use new KL library & new prior names, and compute KLs for the 5 new latent random variables