xFit measures what a difference image holds at each candidate. Its nonlinear source and dipole fits use a batched Levenberg-Marquardt least-squares solver; fixed-position photometry uses a direct linear solution. A dipole is a positive and a negative lobe side by side. A source that moved between the two images leaves one, and so does a slightly misaligned star.
You need: cuPhoton installed with the CPU profile (Chapter 1), and a difference image with its variance and candidate positions (Chapter 6). This page makes them from a sample pair.
NoteSample data
This page runs on a small synthetic pair. The printed values only demonstrate the API.
Create the sample pair (NumPy and Astropy only)
# Writes sample/reference.fits and sample/target.fits: two synthetic star fields# of the same sky on different pixel grids, with different blur and sky level.# The target also holds one new source and one star three times brighter.# Noise is independent Gaussian noise with variance expected_image / gain;# this teaching sample has no read noise, saturation, or detector artifacts.from pathlib import Pathimport numpy as npfrom astropy.io import fitsfrom astropy.wcs import WCSrng = np.random.default_rng(7)def make_wcs(shape, scale_arcsec, rotation_deg, crpix_shift): w = WCS(naxis=2) w.wcs.ctype = ["RA---TAN", "DEC--TAN"] w.wcs.crval = [150.0, 2.0] w.wcs.crpix = [shape[1] /2+ crpix_shift[0], shape[0] /2+ crpix_shift[1]] c, s = np.cos(np.radians(rotation_deg)), np.sin(np.radians(rotation_deg)) w.wcs.cd = scale_arcsec /3600* np.array([[-c, s], [s, c]])return wdef render(shape, wcs, ra, dec, flux, sigma, sky, gain): x, y = wcs.world_to_pixel_values(ra, dec) yy, xx = np.mgrid[: shape[0], : shape[1]] image = np.full(shape, sky, dtype=float)for xi, yi, fi inzip(x, y, flux): image += fi / (2* np.pi * sigma**2) * np.exp(-((xx - xi) **2+ (yy - yi) **2) / (2* sigma**2))return rng.normal(image, np.sqrt(image / gain)).astype(np.float32)ra =150.0+ rng.uniform(-0.016, 0.016, 150)dec =2.0+ rng.uniform(-0.016, 0.016, 150)flux =10** rng.uniform(2.5, 4.5, 150)target_flux = np.r_[3* flux[0], flux[1:], 4000.0] # star 0 brightens; one new sourcetarget_ra, target_dec = np.r_[ra, 150.004], np.r_[dec, 2.003]Path("sample").mkdir(exist_ok=True)for name, shape, wcs, sky_coords, fl, sigma, sky, gain in [ ("reference", (300, 300), make_wcs((300, 300), 0.26, 0.0, (0, 0)), (ra, dec), flux, 1.2, 100.0, 8.0), ("target", (256, 256), make_wcs((256, 256), 0.27, 3.0, (4.3, -2.7)), (target_ra, target_dec), target_flux, 1.8, 120.0, 2.0),]: header = wcs.to_header() header["GAIN"] = (gain, "electrons per data number (DN)") header["PSFSIG"] = (sigma, "Gaussian PSF sigma [px]") fits.writeto(f"sample/{name}.fits", render(shape, wcs, *sky_coords, fl, sigma, sky, gain), header, overwrite=True)
7.1 Make a difference image
The two blocks run the alignment and subtraction of Chapter 6 in short form, without the refit.
import numpy as npfrom astropy.wcs import WCSfrom cuphoton.core.fits_io import read_fits_imagesfrom cuphoton.xrep import BBox, ReprojectionSpec, make_wcs_pair_transform, reproject_masked_arraytarget_read, reference_read = (read_fits_images(f"sample/{name}.fits", [0]) for name in ("target", "reference"))target, reference = (r.arrays[0].astype(np.float64) for r in (target_read, reference_read))target_header, reference_header = target_read.infos[0].header, reference_read.infos[0].headertarget_variance = np.clip(target, 0, None) / target_header["GAIN"]reference_variance = np.clip(reference, 0, None) / reference_header["GAIN"]spec = ReprojectionSpec(mapping=make_wcs_pair_transform(WCS(target_header), WCS(reference_header)), output_bbox=BBox(0, 0, target.shape[1], target.shape[0])) # (x, y) origin, then the sizealigned = reproject_masked_array(reference, reference_variance, np.zeros(reference.shape, dtype=bool), spec, backend="cpu")
The two candidates are the brightened star and the new source, at (y, x) pixels.
7.2 Cut cutouts around the candidates
xFit takes cutouts as an array of shape (batch, y, x). window_stamps from xPois cuts a square window around each (x, y) position, with NaN outside the image. Cut the difference and its variance the same way.
from cuphoton.xpois import window_stampspositions = np.array([region["peak_yx"][::-1] for region in regions], dtype=float) # (x, y)cutouts = window_stamps(difference, positions, 31)cutout_variance = window_stamps(difference_variance, positions, 31)usable = np.isfinite(cutouts) & np.isfinite(cutout_variance) & (cutout_variance >0)print(cutouts.shape, usable.sum(axis=(1, 2)), "usable pixels")
(2, 31, 31) [868 961] usable pixels
The first cutout reaches into the edge band and has fewer usable pixels. Inside a cutout, xFit positions are (x, y) pixels from the center pixel, which is (0, 0). The cutouts are 31 pixels wide to leave room for the quality flags’ optional noise calibration (Section 7.9); the flags themselves also work on smaller cutouts.
7.3 Choose the call
xFit has four ways to measure a cutout:
Call
Fits
Returns
Use when
fit_single_sources
one Gaussian lobe: position, amplitude and, unless fixed, widths and orientation
SingleSourceFitResult
a source is in one image only: new, brightened or faded
forced_psf_photometry
the flux of a PSF at a fixed position
ForcedPhotometryResult
you know the position and the PSF, and want the flux
fit_dipoles with GaussianDipoleModel
two Gaussian lobes of opposite sign with one amplitude and shape
DipoleFitResult
a candidate can hold two lobes, and the PSF is close to Gaussian
fit_dipoles with StampDipoleModel
two shifted copies of a sampled PSF image with one flux
DipoleFitResult
you have a PSF image, and the PSF is not Gaussian
derived_dipole_quantities and dipole_quality_flags read a Gaussian fit_dipoles result.
7.4 Fit one source per cutout
fit_single_sources fits one Gaussian lobe per cutout, from the start position you pass. With psf_sigma the lobe is round with that width, and only its amplitude and position are fitted. Without it, the widths and orientation are free. sign=+1 starts a brightening source and sign=-1 a fading one; without sign, a matched filter at the start sets it. This initializes the amplitude, whose fitted sign is free to change.
from cuphoton.xfit import fit_single_sourcespsf_sigma = target_header["PSFSIG"] # the target PSF width, in pixelsstart = np.zeros((len(positions), 2)) # (x, y) offsets from each cutout center, all zerosources = fit_single_sources(cutouts, cutout_variance, start, psf_sigma, mask=usable, sign=1)for i inrange(len(sources.x)):print(f"x {sources.x[i]:+.2f} +- {sources.position_error_x[i]:.2f}, y {sources.y[i]:+.2f} +- "f"{sources.position_error_y[i]:.2f}, flux {sources.flux[i]:.0f} +- {sources.flux_error[i]:.0f}, "f"converged {sources.converged[i]}")
x +0.46 +- 0.02, y -0.02 +- 0.02, flux 12234 +- 114, converged True
x +0.20 +- 0.04, y +0.05 +- 0.04, flux 4140 +- 71, converged True
flux is the nonnegative lobe-flux magnitude, \(2\pi|A|\sigma_x\sigma_y\) for fitted amplitude \(A\). Recover a signed difference flux with np.sign(sources.parameters[:, 0]) * sources.flux; sources.sign records the starting sign. The result also holds parameters, covariance, standard_errors, chi_square, uncertainty_valid and residuals (model minus data). Its covariance is a local, first-order estimate assuming independent pixels and the supplied variances. It does not account for a wrong PSF or pixel correlations.
7.5 Measure the flux at a fixed position
forced_psf_photometry fits only the flux of a unit-flux Gaussian PSF at fixed positions. The model is linear in the flux. The flux and its error therefore follow without iteration. Pass the fitted positions, or positions from a catalog.
from cuphoton.xfit import forced_psf_photometryfixed = np.column_stack([sources.x, sources.y]) # (x, y) from each cutout centerphotometry = forced_psf_photometry(cutouts, cutout_variance, fixed, psf_sigma, mask=usable)print("flux", photometry.flux.round(0), "error", photometry.flux_error.round(0), "valid", photometry.valid)
Here flux is signed: a fading source gives a negative value. flux_error assumes that the position and the PSF are right and the pixel noise is independent. It leaves out the uncertainty of a position fitted on the same image. To replace the Gaussian, pass unit-flux PSF cutouts, already evaluated at the positions, as psf=; the function does not normalize them. cuphoton.xpois.forced_flux_at applies this estimator to a full image at (x, y) positions. Set local_sky=False to match this example: its default subtracts each window’s median, without adding the uncertainty of that sky estimate to the flux error.
7.6 Build a dipole cutout
The sample pair has no dipole. Build one with NumPy. A star of flux 6,000 sits at (x, y) = (2, 1) in the target and at (-1, -1) in the reference. The difference is target minus matched reference. The star therefore leaves a positive lobe where it is in the target and a negative lobe where it was in the reference. The noise matches the sky of the difference image.
GaussianDipoleModel draws one elliptical Gaussian twice, with opposite signs. Both lobes share one amplitude, their peak height, and one shape, set by the widths sigma_x and sigma_y and the orientation theta. The other four parameters are the lobe centers x_pos, y_pos, x_neg and y_neg. The lobe at (x_pos, y_pos) has the sign of amplitude, and the lobe at (x_neg, y_neg) has the opposite sign. GaussianDipoleModel((31, 31)) takes the cutout shape, (height, width), which must match the images. model="gaussian" builds the same model from the image shape. Fit the two candidates and the dipole in one batch, with a mask that excludes every non-finite pixel.
Check converged and uncertainty_valid before you use a fit. The parameters printed are the dipole’s. Without initial, the fit starts its lobes at the brightest and faintest pixels. A round PSF leaves theta poorly constrained. fit_dipoles does not order the widths, because (sigma_y, sigma_x, theta + pi/2) describes the same ellipse. The result holds these fields:
Field
Holds
parameters, parameter_names
one row of fitted parameters per cutout
standard_errors, covariance
the parameter errors and their full covariance
status, converged, evaluations
how the solver stopped, and after how many evaluations
the chi-square of a zero model, and the fit’s improvement on it
uncertainty_valid, uncertainty_reason
whether the covariance is usable, and why not
residuals
model minus data, in the cutout shape
backend, device, dtype, solver_version
where and in which precision the fit ran, and the solver rules it used
7.8 Derive catalog quantities
derived_dipole_quantities turns each Gaussian dipole fit into its center, the position of the positive lobe (where the source is in the target), the lobe separation and the lobe flux. It propagates the full covariance to each error. The cross terms matter, because amplitude and widths are anticorrelated.
The errors are NaN where uncertainty_valid is false. They are first order, and apply while they are small next to the widths and the separation.
7.9 Screen the fits with quality flags
A converged fit can still be wrong. dipole_quality_flags returns one integer of DipoleQualityFlag bits per fit. Set width_range from the PSF, for example 0.5 to 2.5 times its width.
from cuphoton.xfit import DipoleQualityFlag, dipole_quality_flagsflags = dipole_quality_flags(dipoles, images, variance=variance, mask=mask, width_range=(0.5* psf_sigma, 2.5* psf_sigma))for i, value inenumerate(flags):print(i, int(value), [flag.name for flag in DipoleQualityFlag if value & flag])
Both candidates set NEGATIVE_LOBE_UNSUPPORTED and LOBE_OUTSIDE_STAMP. They are brightening single sources, with the spare lobe placed outside the cutout. Measure them with fit_single_sources. The new source’s spare lobe ran too far out for the lobe test, which then sets both NEGATIVE_LOBE_UNSUPPORTED and POSITIVE_LOBE_UNSUPPORTED. The dipole sets no flag. Read the bits together:
Bit
Flag
Set when
1
UNCERTAINTY_INVALID
the fit has no valid covariance
2
WIDTH_OUT_OF_RANGE
a fitted width lies outside width_range
4
LOBE_OUTSIDE_STAMP
a lobe center lies outside the cutout
8
NO_LOBE_AT_CANDIDATE
no lobe lies within center_radius pixels of the cutout center
16
LOBES_OVERLAP
the lobes are closer than overlap_separation pixels
32
NEGATIVE_LOBE_UNSUPPORTED
the lobes do not overlap and the negative lobe fails the lobe_snr significance test
64
POSITIVE_LOBE_UNSUPPORTED
the lobes do not overlap and the positive lobe fails the lobe_snr significance test
Test one bit with (flags & int(DipoleQualityFlag.LOBES_OVERLAP)) != 0. Exactly one unsupported-lobe bit suggests a brightening or fading single source; read it with the other flags, since a failed fit or unusable lobe can also set it. Closely overlapping lobes make amplitude and separation trade against each other. With calibrate_noise=True, the default, the flags can reduce lobe significance using its robust spread at empty positions of the cutout. This estimates the effect of correlated noise on the lobe test, not on the fitted covariance. Calibration needs at least 12 usable empty positions; even a large cutout can have too few.
7.10 Look at the result
Figure 7.1 shows the three cutouts with the fitted positions.
Figure 7.1: The three cutouts, each on its own red-blue scale centered on zero, red positive. Left and middle: the brightened star and the new source from the sample difference image, with the fit_single_sources position (black cross). The gray band at the bottom of the left cutout lies outside the valid pixels. Right: the synthetic dipole, with the fitted positive (plus) and negative (circle) lobe centers. Synthetic data, for demonstration only.
7.11 Use a sampled PSF
StampDipoleModel shifts copies of a PSF image instead of drawing Gaussians. Its five parameters are the two lobe centers and one shared flux. Use it when you have an image of the PSF, such as a survey’s PSF model, and the PSF is not Gaussian. The block uses a sampled Gaussian as the PSF image.
At the native pixel size, the default bilinear shift biases sub-pixel offsets. For sub-pixel work, pass a PSF image oversampled by a factor n with scale=1 / n. evaluation also takes "bilinear-vignetted" or "finite-volume". Sampled-PSF fits always use finite differences and run on the numpy and cupy backends only.
7.12 Options
The options of fit_dipoles:
Option
Default
Effect
model
required
"gaussian", a GaussianDipoleModel or a StampDipoleModel
initial
None
starting parameters, one row per cutout. None starts at the brightest and faintest pixels.
mask
None
the pixels to fit. It must exclude non-finite image values and non-finite or nonpositive variances.
variance
None
per-pixel variance. Without it, the covariance is scaled by the reduced chi-square.
mode
"difference"
"split" takes (batch, 3, y, x) cutouts: difference, positive and negative planes, treated as independent
backend
"auto"
numpy, cupy, cutile or native. auto picks CuPy when a GPU is usable, else NumPy.
config
None
an LMConfig with solver tolerances, the evaluation budget and the solver version
fusion
False
fuses the Gaussian model and its derivatives on cupy or cutile
If the split-mode difference is computed from the other two planes, those inputs are correlated. The split-mode degrees of freedom and uncertainties do not account for that dependence. Use difference mode for the statistical interpretation shown here.
The LMConfig fields you change first:
Option
Default
Effect
max_evaluations
None
the evaluation budget of each fit. None allows 200 * (n + 1) for n model parameters, 1,800 for a Gaussian dipole. The flag is --max-evaluations.
solver_version
"v0.2"
the stopping, rank and covariance rules. "v0.1.3" selects the earlier rules for comparisons with earlier results. The flag is --solver-version.
The options of dipole_quality_flags:
Option
Default
Effect
width_range
(0.5, 5.0)
allowed lobe widths, in pixels
center_radius
5.0
how far from the cutout center a lobe can sit, in pixels
overlap_separation
2.0
the lobe distance below which LOBES_OVERLAP is set, in pixels
lobe_snr
3.0
the significance a lobe needs, in standard errors
calibrate_noise
True
rescales lobe significance by the spread at empty positions
7.13 On a GPU
With the gpu profile, pass backend="cupy" to fit_dipoles, which still returns NumPy arrays. fit_dipoles_device keeps the result on the GPU. The other xFit calls return CuPy arrays for CuPy inputs. The cutile and native backends fit Gaussian models only, as the xFit reference describes. cutile needs the cutile extra, and native needs the native extension (Chapter 1).
CuPy flags do not combine with DipoleQualityFlag members. Convert a member with int() before you test a bit.
7.14 Command line
cuphoton xfit fit-dipoles fits a batch of cutouts saved as an NPZ archive. The archive needs candidate_id and images, and can hold initial, mask, variance and stamp_basis.
data-validate prints the array shapes and "valid": true. fit-dipoles writes summary.json, fits.parquet and fit-arrays.npz into a new directory and prints the summary. Its metrics block reads:
fits.parquet holds one row per candidate, with each parameter and its standard error:
import pandas as pdtable = pd.read_parquet("runs/dipoles/fits.parquet")print(table[["candidate_id", "uncertainty_valid", "reduced_chi_square", "x_pos", "y_pos"]].round(2).to_string(index=False))
candidate_id uncertainty_valid reduced_chi_square x_pos y_pos
star True 0.92 0.46 -0.02
new False 0.98 0.20 0.05
dipole True 1.01 2.03 1.05
The command line fits dipoles only. To cut the cutouts from aligned FITS images, pass a JSON candidate manifest as --input. The xFit reference describes its format.
7.15 Troubleshooting
ValueError: images must be finite on included pixels
Pass mask=np.isfinite(images) & np.isfinite(variance) & (variance > 0) so invalid image and variance pixels are both excluded.
UserWarning: 21x21 stamps hold fewer than 12 calibration positions
The optional noise calibration needs at least 12 usable empty positions. Square cutouts of 26 pixels or less cannot provide them. Cut larger cutouts, or pass calibrate_noise=False and interpret the lobe tests without a correlation correction.
uncertainty_valid is False
Read uncertainty_reason. A rank-deficient Jacobian means the data do not constrain a parameter, such as a spare lobe on a single source. Check the quality flags.
status is max_evaluations
The fit used up its evaluation budget before it converged. A slowly converging fit can need many evaluations under the default "v0.2" rules. converged is then False and the uncertainties are invalid. Raise max_evaluations in LMConfig, or pass --max-evaluations. A larger budget takes longer and does not guarantee convergence.
ValueError: derived quantities require a Gaussian dipole fit
derived_dipole_quantities and dipole_quality_flags take GaussianDipoleModel fits only.
ValueError: fusion requires backend='cupy' or 'cutile'
Fusion runs on cupy or cutile, with the Gaussian model only.
ERROR: output directory already exists
fit-dipoles needs a new --output-dir for every run.
