commit - c42a195b1cde4822b9b0c4f2f86d627a2c9cac15
commit + 6360d2173215c81e8cf12c81a99bb96eb23dce9b
blob - 60743bf6ee3c9ac94d77ef2fd17cb30464776461
blob + aa698b9bafa0e2485b4dfbb51ecb14bbaf776317
--- Project.toml
+++ Project.toml
[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"
[compat]
CSV = "0.10.16"
+Distributed = "1.11.0"
FFTW = "1.10.0"
FITSIO = "0.17.5"
HTTP = "2.6.4"
blob - 94f982af48c65fd43a71f7b22f2d6e59d0aed721
blob + 73597de260b27d8f7525c16d71419d0145c78157
--- docs/src/design-refinements.md
+++ docs/src/design-refinements.md
`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 on the same real 30-frame field-451
+reference set, 8 worker processes: no crash, and a real 331.8s vs. the
+~720s sequential baseline above — about **2.2x**, 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
using Interpolations
using FFTW
using LsqFit
+using Distributed
include("detection.jl")
include("linking.jl")
blob - 57cc7c697c1903ae896981cd5f18e42f911f8aaf
blob + 61d7ae6f5b300a3eff0b37be6fc96410ccd97205
--- src/reference.jl
+++ src/reference.jl
"""
- 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
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 on real data (real 30-frame ZTF
+field-451 reference set, 8 worker processes) to run without crashing,
+in 331.8s vs. an ~720s sequential baseline (the ~24s/frame figure above)
+— about 2.2x, 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
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)
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)
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
using Test
using Statistics
using Reproject
+using Distributed
@testset "AsteroidPipeline.jl" begin
@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
