diff --git a/benchmark/accuracy/run_celeste_on_field.jl b/benchmark/accuracy/run_celeste_on_field.jl index e239f840..1524afbc 100644 --- a/benchmark/accuracy/run_celeste_on_field.jl +++ b/benchmark/accuracy/run_celeste_on_field.jl @@ -4,6 +4,7 @@ import ArgParse using DataFrames import Celeste.AccuracyBenchmark +import Celeste.Configs import Celeste.Infer import Celeste.ParallelRun import Celeste.SDSSIO @@ -60,6 +61,7 @@ end neighbor_map = Infer.find_neighbors(target_sources, catalog_entries, images) results = AccuracyBenchmark.run_celeste( + Configs.Config(), catalog_entries, target_sources, images, diff --git a/src/AccuracyBenchmark.jl b/src/AccuracyBenchmark.jl index e4eef2be..2b230f35 100644 --- a/src/AccuracyBenchmark.jl +++ b/src/AccuracyBenchmark.jl @@ -6,6 +6,7 @@ import FITSIO import StaticArrays import WCS +import Celeste.Configs import Celeste.DeterministicVI import Celeste.DeterministicVIImagePSF import Celeste.Infer @@ -743,11 +744,13 @@ end # Run Celeste with any combination of single/joint inference and MOG/FFT model function run_celeste( - catalog_entries, target_sources, images; use_joint_inference=false, use_fft=false + config::Configs.Config, catalog_entries, target_sources, images; + use_joint_inference=false, use_fft=false ) neighbor_map = Infer.find_neighbors(target_sources, catalog_entries, images) if use_joint_inference ParallelRun.one_node_joint_infer( + config, catalog_entries, target_sources, neighbor_map, @@ -761,6 +764,7 @@ function run_celeste( infer_source_callback = DeterministicVI.infer_source end ParallelRun.one_node_single_infer( + config, catalog_entries, target_sources, neighbor_map, diff --git a/src/Celeste.jl b/src/Celeste.jl index 5327c43b..cadcabcb 100644 --- a/src/Celeste.jl +++ b/src/Celeste.jl @@ -3,6 +3,7 @@ __precompile__() module Celeste # submodules +include("Configs.jl") include("Log.jl") include("SensitiveFloats.jl") diff --git a/src/Configs.jl b/src/Configs.jl new file mode 100644 index 00000000..a45cb664 --- /dev/null +++ b/src/Configs.jl @@ -0,0 +1,14 @@ +module Configs + +type Config + # A minimum pixel radius to be included around each source. + min_radius_pix::Float64 + + function Config() + config = new() + config.min_radius_pix = 8.0 + config + end +end + +end diff --git a/src/DeterministicVI.jl b/src/DeterministicVI.jl index c88f03d7..ccdecc2b 100644 --- a/src/DeterministicVI.jl +++ b/src/DeterministicVI.jl @@ -4,6 +4,8 @@ Calculate value, gradient, and hessian of the variational ELBO. module DeterministicVI using Base.Threads: threadid, nthreads + +import ..Configs using ..Model import ..Model: BivariateNormalDerivatives, BvnComponent, GalaxyCacheComponent, GalaxySigmaDerivs, SkyPatch, @@ -109,10 +111,10 @@ Arguments: neighbors: the other light sources near `entry` entry: the source to infer """ -function infer_source(images::Vector{Image}, +function infer_source(config::Configs.Config, + images::Vector{Image}, neighbors::Vector{CatalogEntry}, - entry::CatalogEntry; - min_radius_pix=Nullable{Float64}()) + entry::CatalogEntry) if length(neighbors) > 100 msg = string("objid $(entry.objid) [ra: $(entry.pos)] has an excessive", "number ($(length(neighbors))) of neighbors") @@ -125,11 +127,18 @@ function infer_source(images::Vector{Image}, cat_local = vcat([entry], neighbors) vp = init_sources([1], cat_local) patches = Infer.get_sky_patches(images, cat_local) - Infer.load_active_pixels!(images, patches, min_radius_pix=min_radius_pix) + Infer.load_active_pixels!(config, images, patches) ea = ElboArgs(images, vp, patches, [1]) f_evals, max_f, max_x, nm_result = NewtonMaximize.maximize!(elbo, ea) return vp[1] end +# legacy wrapper +function infer_source(images::Vector{Image}, + neighbors::Vector{CatalogEntry}, + entry::CatalogEntry) + infer_source(Configs.Config(), images, neighbors, entry) +end + end diff --git a/src/DeterministicVIImagePSF.jl b/src/DeterministicVIImagePSF.jl index 437fa344..770a61fc 100644 --- a/src/DeterministicVIImagePSF.jl +++ b/src/DeterministicVIImagePSF.jl @@ -7,6 +7,8 @@ module DeterministicVIImagePSF using StaticArrays, DiffBase, Compat +import ..Configs + import ..DeterministicVI: ElboArgs, ElboIntermediateVariables, StarPosParams, GalaxyPosParams, CanonicalParams, VariationalParams, @@ -46,14 +48,14 @@ export elbo_likelihood_with_fft!, FSMSensitiveFloatMatrices, FFTElboFunction, load_fsm_mat -function infer_source_fft(images::Vector{Image}, +function infer_source_fft(config::Configs.Config, + images::Vector{Image}, neighbors::Vector{CatalogEntry}, - entry::CatalogEntry; - min_radius_pix=Nullable{Float64}()) + entry::CatalogEntry) cat_local = vcat([entry], neighbors) vp = init_sources([1], cat_local) patches = get_sky_patches(images, cat_local) - load_active_pixels!(images, patches, min_radius_pix=min_radius_pix) + load_active_pixels!(config, images, patches) ea_fft, fsm_mat = initialize_fft_elbo_parameters(images, vp, patches, [1], use_raw_psf=true) elbo_fft_opt = FFTElboFunction(fsm_mat) @@ -63,6 +65,12 @@ function infer_source_fft(images::Vector{Image}, return vp[1] end +# legacy wrapper +function infer_source_fft(images::Vector{Image}, + neighbors::Vector{CatalogEntry}, + entry::CatalogEntry) + infer_source_fft(Configs.Config(), images, neighbors, entry) +end function infer_source_fft_two_step(images::Vector{Image}, neighbors::Vector{CatalogEntry}, diff --git a/src/GalsimBenchmark.jl b/src/GalsimBenchmark.jl index 4d433b57..16756252 100644 --- a/src/GalsimBenchmark.jl +++ b/src/GalsimBenchmark.jl @@ -5,11 +5,13 @@ module GalsimBenchmark using DataFrames import FITSIO -import Celeste: AccuracyBenchmark, Infer +import Celeste.AccuracyBenchmark +import Celeste.Configs +import Celeste.Infer const GALSIM_BENCHMARK_DIR = joinpath(Pkg.dir("Celeste"), "benchmark", "galsim") const LATEST_FITS_FILENAME_DIR = joinpath(GALSIM_BENCHMARK_DIR, "latest_filenames") -const ACTIVE_PIXELS_MIN_RADIUS_PX = Nullable(40.0) +const ACTIVE_PIXELS_MIN_RADIUS_PX = 40.0 function get_latest_fits_filename(label) latest_fits_filename_holder = joinpath( @@ -66,10 +68,6 @@ end function run_benchmarks(; test_case_names=String[], print_fn=println, joint_inference=false, use_fft=false) - function infer_source_min_radius(args...; kwargs...) - infer_source_callback(args...; min_radius_pix=ACTIVE_PIXELS_MIN_RADIUS_PX, kwargs...) - end - latest_fits_filename = get_latest_fits_filename("galsim_benchmarks") full_fits_path = joinpath(GALSIM_BENCHMARK_DIR, "output", latest_fits_filename) extensions = AccuracyBenchmark.read_fits(full_fits_path) @@ -89,7 +87,10 @@ function run_benchmarks(; test_case_names=String[], print_fn=println, joint_infe truth_catalog_df = extract_catalog_from_header(header) catalog_entries = AccuracyBenchmark.make_initialization_catalog(truth_catalog_df, false) target_sources = collect(1:num_sources) + config = Configs.Config() + config.min_radius_pix = ACTIVE_PIXELS_MIN_RADIUS_PX results = AccuracyBenchmark.run_celeste( + config, catalog_entries, target_sources, images, diff --git a/src/Infer.jl b/src/Infer.jl index 8a64b184..93a6f66e 100644 --- a/src/Infer.jl +++ b/src/Infer.jl @@ -12,6 +12,7 @@ module Infer import WCS using StaticArrays +import ..Configs using ..Model import ..Log @@ -127,14 +128,13 @@ objective function. Non-standard arguments: noise_fraction: The proportion of the noise below which we will remove pixels. - min_radius_pix: A minimum pixel radius to be included. """ -function load_active_pixels!(images::Vector{Image}, +function load_active_pixels!(config::Configs.Config, + images::Vector{Image}, patches::Matrix{SkyPatch}; exclude_nan=true, - noise_fraction=0.5, - min_radius_pix=Nullable{Float64}()) - min_radius_pix = get(min_radius_pix, 8.0) + noise_fraction=0.5) + @show config.min_radius_pix S, N = size(patches) for n = 1:N, s=1:S @@ -155,7 +155,7 @@ function load_active_pixels!(images::Vector{Image}, # include pixels that are close, even if they aren't bright sq_dist = (h - p.pixel_center[1])^2 + (w - p.pixel_center[2])^2 - if sq_dist < min_radius_pix^2 + if sq_dist < config.min_radius_pix^2 p.active_pixel_bitmap[h2, w2] = true continue end @@ -173,6 +173,19 @@ function load_active_pixels!(images::Vector{Image}, end end +# legacy wrapper +function load_active_pixels!(images::Vector{Image}, + patches::Matrix{SkyPatch}; + exclude_nan=true, + noise_fraction=0.5) + load_active_pixels!( + Configs.Config(), + images, + patches, + exclude_nan=exclude_nan, + noise_fraction=noise_fraction, + ) +end # The range of image pixels in a vector of patches function get_active_pixel_range( diff --git a/src/ParallelRun.jl b/src/ParallelRun.jl index 9487a7a8..bf02ef3c 100644 --- a/src/ParallelRun.jl +++ b/src/ParallelRun.jl @@ -5,6 +5,7 @@ using Base.Threads import FITSIO import JLD +import ..Configs import ..Log using ..Model import ..SDSSIO @@ -186,7 +187,8 @@ end """ Optimize the `ts`th element of `sources`. """ -function process_source(ts::Int, +function process_source(config::Configs.Config, + ts::Int, target_sources::Vector{Int}, catalog::Vector{CatalogEntry}, neighbor_map::Vector{Vector{Int}}, @@ -198,7 +200,7 @@ function process_source(ts::Int, neighbors = catalog[neighbor_map[ts]] tic() - vs_opt = infer_source_callback(images, neighbors, entry) + vs_opt = infer_source_callback(config, images, neighbors, entry) ntputs(nodeid, threadid(), "processed objid $(entry.objid) in $(toq()) secs") return OptimizedSource(entry.thing_id, entry.objid, @@ -212,7 +214,8 @@ end Use multiple threads to process each target source with the specified callback and write the results to a file. """ -function one_node_single_infer(catalog::Vector{CatalogEntry}, +function one_node_single_infer(config::Configs.Config, + catalog::Vector{CatalogEntry}, target_sources::Vector{Int}, neighbor_map::Vector{Vector{Int}}, images::Vector{Image}; @@ -238,7 +241,7 @@ function one_node_single_infer(catalog::Vector{CatalogEntry}, end try - result = process_source(ts, target_sources, catalog, neighbor_map, + result = process_source(config, ts, target_sources, catalog, neighbor_map, images; infer_source_callback=infer_source_callback) @@ -268,6 +271,23 @@ function one_node_single_infer(catalog::Vector{CatalogEntry}, return results end +# legacy wrapper +function one_node_single_infer(catalog::Vector{CatalogEntry}, + target_sources::Vector{Int}, + neighbor_map::Vector{Vector{Int}}, + images::Vector{Image}; + infer_source_callback=infer_source, + timing=InferTiming()) + one_node_single_infer( + Configs.Config(), + catalog, + target_sources, + neighbor_map, + images, + infer_source_callback=infer_source_callback, + timing=timing, + ) +end """ Use mulitple threads on one node to fit the Celeste model to sources in a given diff --git a/src/deterministic_vi_image_psf/elbo_image_psf.jl b/src/deterministic_vi_image_psf/elbo_image_psf.jl index 2497a529..c3ef598b 100644 --- a/src/deterministic_vi_image_psf/elbo_image_psf.jl +++ b/src/deterministic_vi_image_psf/elbo_image_psf.jl @@ -391,6 +391,7 @@ end function initialize_fft_elbo_parameters( + config::Configs.Config, images::Vector{Image}, vp::VariationalParams{Float64}, patches::Matrix{SkyPatch}, @@ -398,10 +399,9 @@ function initialize_fft_elbo_parameters( use_raw_psf=true, use_trimmed_psf=true, allocate_fsm_mat=true, - min_radius_pix=Nullable{Float64}()) - +) ea = ElboArgs(images, vp, patches, active_sources, psf_K=1) - load_active_pixels!(images, ea.patches; exclude_nan=false, min_radius_pix=min_radius_pix) + load_active_pixels!(config, images, ea.patches; exclude_nan=false) fsm_mat = nothing if allocate_fsm_mat @@ -411,6 +411,26 @@ function initialize_fft_elbo_parameters( ea, fsm_mat end +# legacy wrapper +function initialize_fft_elbo_parameters( + images::Vector{Image}, + vp::VariationalParams{Float64}, + patches::Matrix{SkyPatch}, + active_sources::Vector{Int}; + use_raw_psf=true, + use_trimmed_psf=true, + allocate_fsm_mat=true, +) + initialize_fft_elbo_parameters( + images, + vp, + patches, + active_sources, + use_raw_psf=use_raw_psf, + use_trimmed_psf=use_trimmed_psf, + allocate_fsm_mat=allocate_fsm_mat, + ) +end @doc """ Return a function callback for an FFT elbo. diff --git a/src/joint_infer.jl b/src/joint_infer.jl index 2bbd3e55..8cfdf7ed 100644 --- a/src/joint_infer.jl +++ b/src/joint_infer.jl @@ -230,14 +230,13 @@ Returns: - Vector of OptimizedSource results """ -function one_node_joint_infer(catalog, target_sources, neighbor_map, images; +function one_node_joint_infer(config::Configs.Config, catalog, target_sources, neighbor_map, images; cyclades_partition::Bool=true, use_fft::Bool=false, batch_size::Int=400, within_batch_shuffling::Bool=true, n_iters::Int=3, - trust_region::NewtonTrustRegion{Float64}=NewtonTrustRegion(), - min_radius_pix::Nullable{Float64}=Nullable{Float64}()) + trust_region::NewtonTrustRegion{Float64}=NewtonTrustRegion()) # Seed random number generator to ensure the same results per run. srand(42) @@ -280,10 +279,9 @@ function one_node_joint_infer(catalog, target_sources, neighbor_map, images; tic() thread_initialize_sources_assignment::Vector{Vector{Vector{Int64}}} = partition_equally(n_threads, n_sources) - initialize_elboargs_sources!(ea_vec, cfg_vec, thread_initialize_sources_assignment, + initialize_elboargs_sources!(config, ea_vec, cfg_vec, thread_initialize_sources_assignment, catalog, target_sources, neighbor_map, images, - trust_region, use_fft, min_radius_pix, - target_source_variational_params) + trust_region, use_fft, target_source_variational_params) Log.info("Done preallocating array of elboargs. Elapsed time: $(toq())") @@ -312,24 +310,47 @@ function one_node_joint_infer(catalog, target_sources, neighbor_map, images; results end -function initialize_elboargs_sources!(ea_vec, cfg_vec, thread_initialize_sources_assignment, +# legacy wrapper +function one_node_joint_infer(catalog, target_sources, neighbor_map, images; + cyclades_partition::Bool=true, + use_fft::Bool=false, + batch_size::Int=400, + within_batch_shuffling::Bool=true, + n_iters::Int=3, + trust_region::NewtonTrustRegion{Float64}=NewtonTrustRegion()) + one_node_joint_infer( + Configs.Config(), + catalog, + target_sources, + neighbor_map, + images, + cyclades_partition=cyclades_partition, + use_fft=use_fft, + batch_size=batch_size, + within_batch_shuffling=within_batch_shuffling, + n_iters=n_iters, + trust_region=trust_region, + ) +end + +function initialize_elboargs_sources!(config::Configs.Config, ea_vec, cfg_vec, + thread_initialize_sources_assignment, catalog, target_sources, neighbor_map, images, - trust_region, use_fft, min_radius_pix, - target_source_variational_params) + trust_region, use_fft, target_source_variational_params) Threads.@threads for i in 1:nthreads() for batch in 1:length(thread_initialize_sources_assignment[i]) - initialize_elboargs_sources_kernel!(ea_vec, cfg_vec, + initialize_elboargs_sources_kernel!(config, ea_vec, cfg_vec, thread_initialize_sources_assignment[i][batch], catalog, target_sources, neighbor_map, images, - trust_region, use_fft, min_radius_pix, + trust_region, use_fft, target_source_variational_params) end end end -function initialize_elboargs_sources_kernel!(ea_vec, cfg_vec, sources, catalog, - target_sources, neighbor_map, images, - trust_region, use_fft, min_radius_pix, +function initialize_elboargs_sources_kernel!(config::Configs.Config, ea_vec, cfg_vec, sources, + catalog, target_sources, neighbor_map, images, + trust_region, use_fft, target_source_variational_params) try for cur_source_index in sources @@ -343,7 +364,7 @@ function initialize_elboargs_sources_kernel!(ea_vec, cfg_vec, sources, catalog, ids_local = vcat([entry_id], neighbor_ids) patches = Infer.get_sky_patches(images, cat_local) - Infer.load_active_pixels!(images, patches, min_radius_pix=min_radius_pix) + Infer.load_active_pixels!(config, images, patches) # Load vp with shared target source params, and also vp # that doesn't share target source params diff --git a/src/multinode_run.jl b/src/multinode_run.jl index aeab543b..6642fc86 100644 --- a/src/multinode_run.jl +++ b/src/multinode_run.jl @@ -1,7 +1,6 @@ using Base.Threads using Gasp - type BoxInfo sources::Vector{Int} catalog::Vector{CatalogEntry} @@ -145,7 +144,14 @@ ntputs(nodeid, tid, "source $ts of $(length(cbox.sources))") # process the source and record the result try - result = process_source(ts, cbox.sources, cbox.catalog, cbox.neighbor_map, cbox.images) + result = process_source( + Configs.Config(), + ts, + cbox.sources, + cbox.catalog, + cbox.neighbor_map, + cbox.images, + ) lock(results_lock) push!(results, result) diff --git a/test/test_eda.jl b/test/test_eda.jl index 25bb7d3a..4f0cce01 100644 --- a/test/test_eda.jl +++ b/test/test_eda.jl @@ -8,8 +8,9 @@ import ..PSF: get_psf_at_point function test_render() images, ea, two_body = SampleData.gen_two_body_dataset(); patches = Infer.get_sky_patches(images, two_body; radius_override_pix=50); - Infer.load_active_pixels!( - images, patches, noise_fraction=Inf, min_radius_pix=Nullable(10)); + config = Configs.Config() + config.min_radius_pix = 10.0 + Infer.load_active_pixels!(config, images, patches, noise_fraction=Inf) ea_fft = DeterministicVIImagePSF.ElboArgs( images, deepcopy(ea.vp), patches, collect(1:ea.S), psf_K=1); n = 3 diff --git a/test/test_galsim_benchmarks.jl b/test/test_galsim_benchmarks.jl index b7f663cd..63885827 100644 --- a/test/test_galsim_benchmarks.jl +++ b/test/test_galsim_benchmarks.jl @@ -5,11 +5,11 @@ using DataFrames import Celeste: GalsimBenchmark GALSIM_CASES_EXERCISED = [ - #"simple_star", + "simple_star", "star_with_noise", "angle_and_axis_ratio_1", "galaxy_with_all", - #"galaxy_with_noise", + "galaxy_with_noise", ] function assert_estimates_are_close(benchmark_results) diff --git a/test/test_infer.jl b/test/test_infer.jl index 602b017e..5890e38c 100644 --- a/test/test_infer.jl +++ b/test/test_infer.jl @@ -1,6 +1,7 @@ ## test the main entry point in Celeste: the `infer` function import JLD +import Celeste.Configs """ test infer with a single (run, camcol, field). @@ -57,7 +58,9 @@ function test_load_active_pixels() # the star is bright but it doesn't cover the whole image. # it's hard to say exactly how many pixels should be active, # but not all of them, and not none of them. - Infer.load_active_pixels!(ea.images, ea.patches; min_radius_pix=Nullable(0.0)) + config = Configs.Config() + config.min_radius_pix = 0 + Infer.load_active_pixels!(config, ea.images, ea.patches) # most star light (>90%) should be recorded by the active pixels num_active_photons = 0.0 @@ -94,7 +97,9 @@ function test_load_active_pixels() # only 2 pixels per image are within 0.6 pixels of the # source's center (10.9, 11.5) - Infer.load_active_pixels!(ea.images, ea.patches; min_radius_pix=Nullable(0.6)) + config = Configs.Config() + config.min_radius_pix = 0.6 + Infer.load_active_pixels!(config, ea.images, ea.patches) for n in 1:ea.N # FIXME: is load active pixels off by (0.5, 0.5)? @@ -106,8 +111,9 @@ end function test_patch_pixel_selection() images, ea, two_body = gen_two_body_dataset(); patches = Infer.get_sky_patches(images, two_body; radius_override_pix=5); - Infer.load_active_pixels!( - images, patches, noise_fraction=Inf, min_radius_pix=Nullable(5)); + config = Configs.Config() + config.min_radius_pix = 5 + Infer.load_active_pixels!(config, images, patches, noise_fraction=Inf) for n in 1:ea.N # Make sure, for testing purposes, that the whole bitmap isn't full.