Tutorial 5

May 19, 2026 · View on GitHub

Introduction

This tutorial performs the same 4-star MCMC fit as Tutorial 4, but using the modern config-driven CLI pipeline introduced in nstarwrap v0.3. Instead of editing Python scripts directly, all parameters live in a single config.toml file. The five CLI commands — nw-prepare, nw-fit, nw-diagnostics, nw-corner, and nw-ds9reg — are always run from the project directory (the one containing config.toml) without any script editing.

What you need in this directory

File	Purpose
`config.toml`	All user parameters
`model.psf`	DAOPHOT PSF model
`4stars.grp`	DAOPHOT group file listing the four stars
`image.fits`	Science image (copy or link or provide path)

The FITS image is not provided here to save space. We will now use the image in: cp ../data_for_tutorials/image.fits.

Conda environments

Make sure to activate the correct environment before starting, and that
pywrapper.so is present (symlink or copy) in this directory — nstarwrap locates it at startup.

Step 0 — Setting up the config.toml configuration file

Everything is controlled by config.toml. Read through it once before running anything. The key sections are explained below.

`[run]`

[run]
run_id = "run0"

Every output file uses this identifier as a prefix. Results go into a subdirectory run0/ that nw-fit creates automatically.

`[data]`

[data]
image_fits = "../data_for_tutorials/image.fits"   # used by nw-diagnostics for display
image_pkl  = "../data_for_tutorials/image.fits.pkl"  # fast binary for DAOPHOT workers
psf_file   = "model.psf"
group_file = "4stars.grp"
box_xmin   = 1070
box_xmax   = 1143
box_ymin   = 1147
box_ymax   = 1213

image_pkl is a fast binary copy of the FITS image (see Step 1 below). box_* are 1-indexed pixel coordinates of the region where DAOPHOT performs the PSF fit.

`[daophot]`

DAOPHOT internal options. fitrad (fitting radius in pixels) controls how far from a star's centroid pixels contribute to the chi-square. These are fixed throughout the MCMC.

`[firstguess]`

[firstguess]
x1 = 1112.536  y1 = 1183.238
x2 = 1112.368  y2 = 1189.931
x3 = 1093.820  y3 = 1171.420
x4 = 1108.780  y4 = 1167.000
flux_fraction1 = 0.657333700306   # f1 = F1 / (F1+F2+F3+F4)
flux_fraction3 = 0.173373372520   # f3 = F3 / (F1+F2+F3+F4)
flux_fraction4 = 0.058225386230   # f4 = F4 / (F1+F2+F3+F4)
err_factor     = 0.122580628

Flux fraction convention. For a 4-star model the 11 free parameters (ndim = 11) are the four (x, y) centroids and three flux fractions. The fourth fraction is always implied: f2 = 1 − f1 − f3 − f4. The code enforces a hard prior that all fractions are positive and their sum is strictly less than 1.

err_factor rescales the pixel variances passed to DAOPHOT. A value different from 1 means the nominal pixel noise is not a perfect description of the data. nw-fit suggests an updated value at the end of the run (look for "Suggested err_factor" in the terminal output).

`[fit]`

The MCMC settings. The most important ones:

Key	Value	Meaning
`ndim`	11	4-star model (8 positions + 3 flux fractions)
`nwalkers`	22	Number of parallel chains; must be ≥ 2 × ndim
`max_n`	10000	Steps per chain
`move_type`	`"stretch"`	emcee StretchMove (same algorithm as Tutorial 4)
`flag_refine_bestfit`	`true`	Polish best MCMC sample with L-BFGS-B after the run
`init_run_id`	`""`	Empty = start walkers from a Gaussian ball around firstguess

`[diagnostics]` and `[corner]`

Visualization parameters. burnin is how many early steps to discard before computing statistics and making plots. thin keeps every n-th sample, reducing memory usage in the plots. For a first look, burnin = 10 and thin = 1 are fine; revisit after inspecting the chain traces (see Step 3).

Step 1 — prepare the image binary (nw-prepare)

nw-prepare converts the FITS image to a fast binary (.pkl) that is memory-mapped efficiently by parallel DAOPHOT workers:

nw-prepare

This creates image.fits.pkl in the current directory.

Step 2 — run the MCMC (nw-fit)

nw-fit

nw-fit reads config.toml, creates the subdirectory run0/, and starts the ensemble sampler. A progress table is printed every chunk steps:

step    accept   log_prob(max)   χ²_min    χ²_r
 1000   0.412    -1263.2         2526.4    1.001
 2000   0.398    -1262.1         2524.2    0.999
  ...

Key output files written to run0/:

File	Contents
`run0.h5`	Full emcee chain (HDFBackend, grows during the run)
`run0.txt`	Terminal log including the "Zeropoint Mag:" line
`run0_nstar.log`	Deduplicated DAOPHOT warnings
`run0_bestfit.json`	Refined best-fit parameters (L-BFGS-B result)
`config.toml`	Snapshot of the exact config used for this run

The .h5 file is written incrementally. You can open a second terminal and run nw-diagnostics or nw-corner while the MCMC is still in progress.

Step 3 — inspect the chains (nw-diagnostics)

nw-diagnostics

Produces a self-contained HTML report in run0/:

run0/run0_diagnostics_001.html
run0/run0_chains_001.png
run0/run0_walkers_001.png

Open the HTML in a browser. The report contains three main panels:

Walker positions. Each walker is drawn as a coloured trajectory on the science image. Circles mark star 1, diamonds mark star 2. If all walkers have converged to tight clusters around each star, mixing is good.

Chain traces. Parameter values as a function of MCMC step, one line per walker. The early transient (burn-in) is visible as a rapid drift before the chains settle into stationary sampling. The point where the chains stabilise is the burn-in length you should set in [diagnostics] and [corner].

Summary table. Best-fit values (★ = refined by L-BFGS-B), 16th/50th/84th percentiles, χ²_min, DOF, χ²_r.

Re-running nw-diagnostics after adjusting burnin writes a new _002.html alongside the first. Old outputs are never overwritten.

Step 4 — corner plot and contour export (nw-corner)

nw-corner

Produces in run0/:

run0/run0_correlation_001.png   ← 11×11 corner plot
run0/run0_summary_001.txt       ← text summary with autocorrelation times
run0/run0_contours_001.json     ← posterior contours for DS9 and nw-compare

The corner plot shows the 2D marginal posterior for every parameter pair. The white circle is the MCMC argmax; the cyan cross is the refined best-fit.

The summary file lists autocorrelation times (τ) for each parameter. Use the trace plots to choose a final burn-in and τ for the thinning (see Tips below).

Step 5 — DS9 region file (nw-ds9reg)

nw-ds9reg

Writes run0/run0.reg. Open it in DS9 alongside the FITS image:

File → Open → image.fits
Region → Load Regions → run0/run0.reg

The region file contains:

White box — the DAOPHOT fitting region.
Coloured polygons — posterior position contours from the contours JSON: red = 1σ, blue = 2σ, magenta = 3σ (one set per star).
Black × markers — best-fit centroids.
Magnitude labels with leader lines — instrumental magnitudes computed from the zeropoint and each star's flux fraction.

Step 6 — compare runs (nw-compare)

Once you have more than one run (e.g., after tuning err_factor or switching move_type to "de" and rerunning with a different run_id), overlay their posteriors in a single HTML page. Update config.toml:

[compare]
input_specs = ["run0", "run1"]
output_stem = "compare"

Then:

nw-compare

The output is a self-contained HTML with interactive σ-level checkboxes in the sidebar and a summary table (χ²_r, DOF, number of samples per run).

For the full list of accepted input_specs formats, run aliases, and figure export options, see the nw-compare section of the pipeline documentation.

Workflow summary

nw-prepare         # once per image (skip here — reusing data_for_tutorials pkl)
     ↓
nw-fit             # MCMC run → run0/
     ↓
nw-diagnostics     # HTML report — check burn-in, update [diagnostics] burnin
nw-corner          # corner plot + contours JSON — check autocorr times
nw-ds9reg          # DS9 region file with posterior contours
nw-compare         # (optional) overlay multiple runs

All commands read config.toml in the current directory. Never edit package source files — only config.toml changes between runs.

Tips

Choosing burn-in

Look at the chain traces in nw-diagnostics. Discard everything before the chains settle into stationary sampling. As a rule of thumb, the burn-in should be at least 2–3 times the autocorrelation time τ printed by nw-corner. Update burnin in both [diagnostics] and [corner] and re-run those commands.