andymiller commented 7 years ago

added support for MCMC within AIS inference

implemented annealed importance sampling to approximate the marginal likelihood of a star vs galaxy
code for running slice-sampling within AIS on both star and galaxy conditional distributions
code for consolidating samples and assembling output dataframes (mcmc/mcmc_infer.jl)
added a process_source_mcmc function in ParallelRun.jl; added a boolean do_vi flag to one_node_single_infer (defaults to VI)

The following call should return a list of res objects (dictionaries)

 res_list = ParallelRun.one_node_single_infer(
     catalog_entries, target_sources, neighbor_map, images;
     config=Config(25.0), do_vi=false)

each object is a dictionary, with the following entries

star_samples and gal_samples : dataframe of parameter samples
star_lls and gal_lls : log joint probability for the star or galaxy samples
ave_pstar : average log probability that the source is a star
star_bootstrap, gal_bootstrap: bootstrap samples of the marginal likelihood of the data under the star model and the gal model, respectively (these are used to compute the posterior prob(star | data)).

The list of dictionaries above can be consolidated into a one-source-per-row dataframe with approximate posterior means and stderr values with the following call

star_summary, gal_summary, joint_summary = MCMC.consolidate_samples(res_list)

Here, star_summary assumes all sources are stars, gal_summary assumes all sources are galaxies, and joint_summary presents star posterior inferences if exp(ave_pstar) > .5, and galaxy posterior inferences otherwise.

Closes #387 . Closes #324 .

codecov[bot] commented 7 years ago

Codecov Report

Merging #716 into master will decrease coverage by 1.2%. The diff coverage is 65.61%.

@@            Coverage Diff             @@
##           master     #716      +/-   ##
==========================================
- Coverage   81.24%   80.03%   -1.21%     
==========================================
  Files          35       37       +2     
  Lines        3642     4033     +391     
==========================================
+ Hits         2959     3228     +269     
- Misses        683      805     +122

Impacted Files	Coverage Δ
src/model/log_prob.jl	`90.57% <ø> (ø)`	:arrow_up:
src/config.jl	`100% <100%> (ø)`	:arrow_up:
src/mcmc/mcmc_infer.jl	`37.83% <37.83%> (ø)`
src/mcmc/mcmc_misc.jl	`48.18% <47.42%> (+14.6%)`	:arrow_up:
src/mcmc/slicesample.jl	`58% <57.89%> (-18.14%)`	:arrow_down:
src/mcmc/mcmc_functions.jl	`69.23% <82.87%> (+23.28%)`	:arrow_up:
src/ParallelRun.jl	`83.01% <85.71%> (+0.09%)`	:arrow_up:
src/mcmc/ais.jl	`92.85% <92.85%> (ø)`
... and 1 more

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 657886d...36c53a7. Read the comment docs.

jeff-regier commented 7 years ago

Thanks for fixing the merge conflicts.

Would you add a unit test that calls ParallelRun.one_node_single_infer with do_vi = false? Maybe use config=Config(2.0) to keep the test fast.

That test should exercise basically all your MCMC code, right? If there's code not exercised by that, can we omit it from the PR?

jeff-regier commented 7 years ago

Stuff like write_star_unit_flux_raw and write_galaxy_unit_flux and render_patch duplicates code in src/Synthetic.jl. We want to avoid duplicating that code especially, since it could lead to us doing slightly different tests.

andymiller commented 7 years ago

I agree there is some overlap in the write_star_unit_flux_raw and render_patch code. I reimplemented the old version from Synthetic to control the offset so I could use a WCS object on a small patch of image. I realize that is now supported in write_star_nmgy! method --- but one thing I'm hesitant of is re-instantiating a SkyPatch (and then loading the active pixels) every time the likelihood function is called.

jeff-regier commented 7 years ago

Would it work to make p (the SkyPatch) an argument to write_star_nmgy!?

andymiller commented 7 years ago

oh also the pixel map in write_star_unit_flux_raw and others is just an array. What if I wrote a general method

function write_star_nmgy!(world_pos::Array{Float64,1}, patch::SkyPatch,  flux::float64, pixels::Matrix{Float64} )
    fs0m = SensitiveFloat{Float64}(length(StarPosParams), 1, false, false)

    H2, W2 = size(p.active_pixel_bitmap)
    for w2 in 1:W2, h2 in 1:H2
        h = p.bitmap_offset[1] + h2
        w = p.bitmap_offset[2] + w2
        Model.star_light_density!(fs0m, p, h, w, pos, false)
        pixels[h, w] += fs0m.v[] * flux
    end
end

and reimplemented what is in Synthetic with

function write_star_nmgy!(img::Image, ce::CatalogEntry)
    p = Model.SkyPatch(img, ce, radius_override_pix=25)
    write_star_nmgy!(ce.pos, p, ce.star_fluxes[img.b], img.pixels) 
end

and did the same for write_galaxy_nmgy!. Then I can call all of this code from the MCMC module.

jeff-regier commented 7 years ago

I just fixed a minor merge conflict with master: please pull before editing.

Looks like make_position_transformations and sample_colors and sample_logr and parameters_from_catalog aren't covered by tests. Are they necessary? Could we safely omit them from the master branch?

No need to write thorough tests, but it's nice if every function in the src directory of master gets called once, just to catch syntax errors if anyone is making changes without totally understanding what's going on.

jeff-regier commented 7 years ago

Ah, this latest merge had more conflicts that I expected. Let's get this PR committed as is, to avoid more merge conflicts.

Andy, please do address the stuff in this thread in another PR. Please start a new branch off of master rather than continuing an old branch---otherwise all the old commits get mixed in.

andymiller commented 7 years ago

ok, i moved all my recent changes to a new branch ais-benchmark and squashed the changes. Will open a new PR shortly.

jeff-regier / Celeste.jl

ais inference #716

Codecov Report