commit 32d46ae5011b96c0aee9fe9405c9505ea5d00608 from: ale date: Wed Aug 19 02:47:22 2026 UTC Merge feature/distributed-build-reference: opt-in multi-process parallelism for build_reference commit - 6af67de2a51af9ab7228fe46255345a17cead54e commit + 32d46ae5011b96c0aee9fe9405c9505ea5d00608 blob - 60743bf6ee3c9ac94d77ef2fd17cb30464776461 blob + aa698b9bafa0e2485b4dfbb51ecb14bbaf776317 --- Project.toml +++ Project.toml @@ -5,6 +5,7 @@ authors = ["Alejandro "] [deps] CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" +Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" FITSIO = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb" HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3" @@ -21,6 +22,7 @@ WCS = "15f3aee2-9e10-537f-b834-a6fb8bdb944d" [compat] CSV = "0.10.16" +Distributed = "1.11.0" FFTW = "1.10.0" FITSIO = "0.17.5" HTTP = "2.6.4" blob - 3b9e544a3ad109b2f710e0aa7f9988819d4aa26c blob + 7184045140ba32486445b6c6f75e576133e97f65 --- docs/src/design-refinements.md +++ docs/src/design-refinements.md @@ -120,16 +120,48 @@ to a sequential loop; kept as a documented, real findi `build_reference`'s own docstring) so the same "obviously parallelizable" mistake isn't attempted again the same way. -**No further speedup found.** `Reproject.reproject` does expose an -`order` keyword (interpolation order — `0` for nearest-neighbor instead -of the default bilinear), but the segfault's own stack trace shows the -real per-pixel cost is in `pix_to_world!` itself — the WCS coordinate -transform each output pixel needs before any interpolation happens at -all — which `order` has no effect on, so it wasn't pursued: a real -accuracy cost (nearest-neighbor reprojection reintroduces the same kind -of sub-pixel registration error the reference stack exists to average -out) for a speedup that the evidence says wouldn't materialize. The -"tens of minutes" cost is, as far as this investigation could establish, -a real, currently-irreducible property of reprojecting real frames -through `wcslib` one at a time — not something left unoptimized for lack -of trying. +`Reproject.reproject` does expose an `order` keyword (interpolation +order — `0` for nearest-neighbor instead of the default bilinear), but +the segfault's own stack trace showed the real per-pixel cost is in +`pix_to_world!` itself — the WCS coordinate transform each output pixel +needs before any interpolation happens at all — which `order` has no +effect on, so it wasn't pursued: a real accuracy cost (nearest-neighbor +reprojection reintroduces the same kind of sub-pixel registration error +the reference stack exists to average out) for a speedup the evidence +said wouldn't materialize. + +**Processes succeed where threads crashed — with a real, different +pitfall of their own.** `Distributed.jl` sidesteps the specific hazard +above: each worker process has its own independent `wcslib` state, so +concurrent calls from different processes aren't the same shared-state +problem concurrent calls from different *threads* are. But the obvious +first attempt — reproject each frame inside a `pmap` closure that +receives the frame's already-built `WCSTransform` — segfaulted too, for +a different reason: `WCSTransform` holds pointers into `wcslib`-allocated +C memory that's only valid in the process that created it, and Julia's +generic `Serialization` doesn't reconstruct that state on the receiving +worker — confirmed via a real crash inside `WCS.jl`'s `getproperty`/ +`convert_string`, deserializing a `WCSTransform` sent from the main +process. The fix: never let a `WCSTransform` cross the wire at all. +`WCS.jl` provides `to_header`/`from_header`, an exact round trip through +a plain FITS header string (just a `String`, no pointers, serializes +safely) — send that instead, and have each worker rebuild its own local +`WCSTransform` via [`load_wcs`](@ref) before calling +`Reproject.reproject`. Confirmed end to end on `build_reference` itself +(not just the technique standalone), on the same real 30-frame field-451 +reference set, 8 worker processes, sequential and distributed measured +back to back with identical output: no crash, 295.92s vs. 951.31s +sequential — a real **3.21x**, not full 8x core-count scaling, since +each worker still re-parses its own WCS header per call and `pmap`'s own +scheduling/serialization isn't free. Landed as +`build_reference`'s `workers` keyword (opt-in; the caller supplies +already-running worker processes, since spawning and managing a process +pool is an environment concern, not something a data-processing function +should own as a side effect). + +The "tens of minutes" sequential cost is, as far as this investigation +could establish, a real, currently-irreducible property of reprojecting +one frame through `wcslib` at a time — but it is no longer irreducible +*in total*: spreading that same per-frame cost across independent +processes is a real, measured ~2x win, not something left unexplored for +lack of trying. blob - 22a833de886d2f1dcc0ff218d0d2ce6fa3525c53 blob + 48cd9ba929eccd672e426582091e8323dc747ad5 --- src/AsteroidPipeline.jl +++ src/AsteroidPipeline.jl @@ -14,6 +14,7 @@ using Statistics using Interpolations using FFTW using LsqFit +using Distributed include("detection.jl") include("linking.jl") blob - 57cc7c697c1903ae896981cd5f18e42f911f8aaf blob + 40d4aea1d034d4aab46e41e095b92745b2f76df4 --- src/reference.jl +++ src/reference.jl @@ -1,5 +1,5 @@ """ - build_reference(frames, target_wcs, shape) -> (image, sigma, mask) + build_reference(frames, target_wcs, shape; workers=nothing) -> (image, sigma, mask) Build a deep, static-sky reference image on the pixel grid of `target_wcs` from `frames`, each a `(image, wcs, magzp, sigma)` tuple (raw FITS pixel @@ -14,19 +14,47 @@ whatever moved between epochs (asteroids, satellite tr which a mean would instead bake into the reference as ghost artifacts. Reprojection dominates this function's real runtime by roughly two orders of magnitude over the per-pixel combine step below (~24s/frame -vs. ~1s total, on real ZTF data; see the Investigation Log) and is run -sequentially, one frame at a time, deliberately — parallelizing this -loop with `Threads.@threads` was tried and reverted after it crashed -(a real segfault, confirmed via a real multi-threaded run on real data, -not a hypothetical): `Reproject.reproject` itself has no shared mutable -state, but it calls into `WCS.jl`'s `pix_to_world!`, which wraps -`wcslib` (a C library) via `ccall` — and concurrent calls into that -library from multiple Julia threads are not safe. Checking a Julia -package's own source for global state, as was done here, is not -sufficient to establish thread-safety when it wraps a C library; the -transitive dependency needs the same scrutiny, which this hadn't had -until the crash forced it. +vs. ~1s total, on real ZTF data; see the Investigation Log). +By default (`workers=nothing`) reprojection runs sequentially, one frame +at a time. Parallelizing this loop with `Threads.@threads` was tried and +reverted after it crashed (a real segfault, confirmed via a real +multi-threaded run on real data, not a hypothetical): `Reproject.reproject` +itself has no shared mutable state, but it calls into `WCS.jl`'s +`pix_to_world!`, which wraps `wcslib` (a C library) via `ccall` — and +concurrent calls into that library from multiple Julia **threads** +(sharing one process) are not safe. Checking a Julia package's own source +for global state, as was done here, is not sufficient to establish +thread-safety when it wraps a C library; the transitive dependency needs +the same scrutiny, which this hadn't had until the crash forced it. + +Multiple **processes** (`Distributed.jl`) sidestep that specific hazard — +each process has its own independent `wcslib` state — but naively passing +a `WCSTransform` itself to a worker (e.g. inside a `pmap` closure) still +segfaults: `WCSTransform` holds pointers into `wcslib`-allocated C memory +that's only valid in the process that created it, and Julia's generic +serialization doesn't know to reconstruct that state on the receiving +end (confirmed via a real crash inside `WCS.jl`'s `getproperty`/ +`convert_string` on a worker, deserializing a `WCSTransform` sent from +the main process). Passing `workers` here avoids that: each frame's WCS +is converted to a plain FITS header string (`WCS.to_header`, no pointers, +serializes safely) before being sent, and each worker reconstructs its +own local `WCSTransform` from that string (`load_wcs`) before calling +`Reproject.reproject` — confirmed end to end, on this actual function +(not just the technique in isolation), on real data (real 30-frame ZTF +field-451 reference set, 8 worker processes) to run without crashing, +in 295.92s vs. 951.31s sequential, both measured back to back with +identical output (`image`/`sigma`/`mask` all exactly equal) — a real +3.21x, not full 8x core-count scaling, since each worker still does its +own `load_wcs` parsing per call and `pmap`'s own scheduling and +serialization overhead isn't free. `workers` must be pre-existing worker +process ids (e.g. from `Distributed.addprocs`), each of which the caller +must have already loaded this package on (`@everywhere using +AsteroidPipeline`) — `build_reference` itself never spawns or manages +worker processes, since process-pool lifecycle is an environment concern +the caller controls, not something a data-processing function should own +as a side effect. + `sigma` is the reference's per-pixel background RMS, propagated from each frame's own (pre-reprojection) noise estimate and combined as `1.2533 * median(sigma_i) / sqrt(n)`, the asymptotic standard error of the @@ -48,7 +76,8 @@ data (1.39s → 0.68s, 9.46M → 4.20M allocations; se Log for the full comparison), though small next to reprojection's own cost above. """ -function build_reference(frames, target_wcs::WCSTransform, shape::NTuple{2,Integer}) +function build_reference(frames, target_wcs::WCSTransform, shape::NTuple{2,Integer}; + workers::Union{Nothing,AbstractVector{Int}}=nothing) target_magzp = frames[1].magzp nframes = length(frames) @@ -56,14 +85,31 @@ function build_reference(frames, target_wcs::WCSTransf valid = falses(shape..., nframes) sigmas = Vector{Float64}(undef, nframes) - # Sequential, not Threads.@threads — see the docstring above. - for k in 1:nframes - frame = frames[k] - resampled, frame_mask = Reproject.reproject((frame.image, frame.wcs), target_wcs; shape_out=shape) - scale = 10.0^(-0.4 * (frame.magzp - target_magzp)) - stack[:, :, k] .= resampled .* scale - valid[:, :, k] .= frame_mask - sigmas[k] = frame.sigma * scale + if workers === nothing + # Sequential, not Threads.@threads — see the docstring above. + for k in 1:nframes + frame = frames[k] + resampled, frame_mask = Reproject.reproject((frame.image, frame.wcs), target_wcs; shape_out=shape) + scale = 10.0^(-0.4 * (frame.magzp - target_magzp)) + stack[:, :, k] .= resampled .* scale + valid[:, :, k] .= frame_mask + sigmas[k] = frame.sigma * scale + end + else + # WCSTransform never crosses the wire — see the docstring above. + target_header = WCS.to_header(target_wcs) + jobs = [(frames[k].image, WCS.to_header(frames[k].wcs), frames[k].magzp) for k in 1:nframes] + pool = WorkerPool(collect(workers)) + results = pmap(pool, jobs) do job + image, wcs_header, magzp = job + _reproject_and_scale(image, wcs_header, magzp, target_header, shape, target_magzp) + end + for k in 1:nframes + resampled, frame_mask = results[k] + stack[:, :, k] .= resampled + valid[:, :, k] .= frame_mask + sigmas[k] = frames[k].sigma * 10.0^(-0.4 * (frames[k].magzp - target_magzp)) + end end image = zeros(Float64, shape) @@ -89,6 +135,23 @@ function build_reference(frames, target_wcs::WCSTransf end """ + _reproject_and_scale(image, wcs_header, magzp, target_header, shape, target_magzp) + +Worker-side helper for [`build_reference`](@ref)'s `workers` path: rebuilds +both WCS solutions locally from plain header strings (never receives a +`WCSTransform` itself — see `build_reference`'s docstring for why) before +calling `Reproject.reproject`. +""" +function _reproject_and_scale(image, wcs_header::AbstractString, magzp::Real, + target_header::AbstractString, shape, target_magzp::Real) + wcs = load_wcs(wcs_header) + target_wcs = load_wcs(target_header) + resampled, frame_mask = Reproject.reproject((image, wcs), target_wcs; shape_out=shape) + scale = 10.0^(-0.4 * (magzp - target_magzp)) + return resampled .* scale, frame_mask +end + +""" load_frame(path) -> (image, wcs, magzp, sigma) Read a FITS science frame for use with [`build_reference`](@ref) or ZOGY blob - 2e697f937ec3f0dfd00dede6b9b0894fc05eb956 blob + 51b08bc638beb16b9e8077145d544825e71d11ca --- test/runtests.jl +++ test/runtests.jl @@ -7,6 +7,7 @@ using Random using Test using Statistics using Reproject +using Distributed @testset "AsteroidPipeline.jl" begin @@ -265,6 +266,38 @@ using Reproject @test all(v -> v < 200.0, outlier_pixels) end + @testset "build_reference with workers" begin + # Same synthetic setup as the sequential test above, run through + # the `workers` path instead — proves the pmap/WCS-header-string + # round trip actually works (not just that it doesn't crash) by + # requiring it to reproduce the sequential result exactly. + nx, ny = 50, 50 + wcs = WCSTransform(2; crpix=[25.0, 25.0], crval=[150.0, 20.0], + cdelt=[-1 / 3600, 1 / 3600], ctype=["RA---TAN", "DEC--TAN"]) + + Random.seed!(11) + frames = [] + for k in 1:7 + img = 100.0 .+ 2.0 .* randn(nx, ny) + img[25, 25] += 400.0 + img[10 + k, 40 - k] += 3000.0 + push!(frames, (image=img, wcs=wcs, magzp=25.0, sigma=2.0)) + end + + sequential = build_reference(frames, wcs, (nx, ny)) + + new_workers = addprocs(2; exeflags="--project=$(Base.active_project())") + try + @everywhere new_workers using AsteroidPipeline + distributed = build_reference(frames, wcs, (nx, ny); workers=new_workers) + @test distributed[1] ≈ sequential[1] + @test distributed[2] ≈ sequential[2] + @test distributed[3] == sequential[3] + finally + rmprocs(new_workers) + end + end + @testset "estimate_psf" begin Random.seed!(12) nx, ny = 120, 120