Imaging (MCMC)

reconstruct_squeeze(data, ft; kwargs...)           -> (image, diagnostics)
reconstruct_squeeze(x_start, data, ft; kwargs...)  -> (image, diagnostics)

SQUEEZE-style MCMC image reconstruction.

The image is represented as a bag of nelements discrete flux quanta on the pixel grid, sampled by Metropolis-Hastings with simulated annealing. This gives positivity and fixed total flux for free, admits non-convex regularizers such as L0, and returns a posterior mean image rather than a MAP point.

x_start selects the starting configuration; omit it for SQUEEZE's own default.

Strings work too ("point_source", "random"). Note the default is a point source, not a random or flat image: the annealing is expected to spread the flux out from it.

ft follows the same convention as reconstruct / reconstruct_bsmem: the image geometry is taken from the Fourier operator, not from keywords.

• Pass a DFT matrix (setup_ft(data, nx, pixsize; mode="dft") or setup_dft) and it is used directly — this is what the sampler wants, and nothing is rebuilt.
• Pass an NFFT plan and nx/pixsize are recovered from it exactly as reconstruct_bsmem does, then a DFT matrix is built internally. The sampler cannot use an NFFT plan: it needs O(nuv) access to a single column per move (the rank-1 update), which a transform-based operator cannot provide.

Keyword arguments

• nelements: number of flux quanta (default default_nelements).
• niter: sweeps; each sweep is nelements proposals.
• nchains: independent annealing restarts, run on separate threads.
• regularizers: OITOOLS-style [["l0", λ], ["tv", λ], ...]. Supported names are l0, tv, entropy, compactness, centering, priorimage. These are SQUEEZE's forms on the integer histogram, not OITOOLS' differentiable versions.
• weights: observable weights, default [1.0,1.0,1.0] = V², T3amp, T3phi, as in reconstruct (padded to 7 internally).
• prior_image: an nx × nx map of per-pixel prior probabilities (C's -p). It becomes an additive penalty -log(p) per quantum, with p <= 0 mapped to 1e12 — which is what makes it a mask: no quantum can ever sit on a zero-prior pixel. Its weight defaults to 1 and can be set explicitly with ["priorimage", λ] in regularizers. Supplying one disables auto-centering, since the prior already fixes the position (C warns about this combination).
• pixsize: only needed to override the value recovered from ft, or to supply one when ft is a bare DFT matrix and outfile is set (the FITS WCS header).
• tmin, chi2_temp, chi2_target: annealing schedule (C's -tm, -ct, -fc).
• f_anywhere, f_copycat: fixed proposal fractions.
• model: an optional SqueezeSparco chromatic star + background model.
• print_every: C-format text diagnostic line every n sweeps.
• monitor: if > 0, redraw a live image and trace display every monitor sweeps (chain 1 only). Off by default, and free when off — one integer comparison per sweep. Enabling it pins chain 1 to the calling thread, because drawing must happen on the main thread; the other chains still run in parallel. monitor_colormap selects the image colour map.
• seed: per-chain streams are a pure function of (seed, chain), so results are reproducible regardless of thread scheduling.

Notes

COPYCAT proposals are asymmetric and irreversible; they are an effective optimizer heuristic for annealing but are not a valid sampler move. Set f_copycat=0 if you want the post-burn-in samples to be a defensible posterior.

reconstruct_squeeze(data, ft; kwargs...)

Start from SQUEEZE's default configuration — a point source at the centre of the grid, the equivalent of running the C code without -i. Identical to passing :point_source as x_start.

reconstruct_squeeze(x_start, data::AbstractArray{<:OIdata}, ft; kwargs...)

Convenience dispatch matching the rest of the reconstruct* family, so the usual data = readoifits(file); ft = setup_ft(data, nx, pixsize; mode="dft") pattern works without indexing. Only the monochromatic, single-epoch case is supported.

reconstruct_squeeze_tempered(data, ft; kwargs...)

Parallel-tempered SQUEEZE reconstruction, built on Pigeons.jl.

Requires Pigeons to be loaded:

using OITOOLS, Pigeons
img, diag = reconstruct_squeeze_tempered(data, ft; n_rounds = 10, n_chains = 10)

Unlike reconstruct_squeeze, which anneals (an optimizer that floors at T = 1), this samples the posterior at β = 1 with a full temperature ladder, and returns a log normalising constant logZ alongside the mean image.

SqueezeSparco(; f_star, ud, env_indx, lambda0, f_bg, bg_indx, free, stepsize)

SPARCO model state for reconstruct_squeeze.

Parameters, in C's order:

free lists the parameters the sampler may vary, by name or index; everything else is held fixed. stepsize is the two-point lattice half-width per parameter and is adapted toward a 0.3 acceptance rate during the run.

The model is only meaningful on polychromatic data: with a single wavelength every (λ/λ0)^x is a constant and the three indices are exactly degenerate.

default_nelements(data, nx)

C's heuristic (squeeze.c:298): 2·ceil(nx · ndata^(1/3)), floored at 500.