Basic Usage

Basic Usage#

How to use fpho with your images.

Getting ready#

First, you will need collect the appropriate imaging data, PSFs, and identify a an initial list or catalog of ‘peaks’ to be fit. A pre-processing script may be used to convert the imaging data into an efficient format useable by forcepho (a pixel-data store and meta-data store). See ./inputs.rst for details

Second, you will need to generate a configuration file with information about data input and output locations and details for the fitting process. See ./configuration.rst for details

A Gaussian mixture approximation to each relevant PSF must be generated, using tools provided with forcepho. These are stored in an HDF5 file, in data groups keyed by FILTER.

Then, the following steps will lead to output that can be post-processed. See ./output.md for details on post-processing

Basic Fitting#

The basic procedure requires several ingredients to be instantiated using the information above. Examples are given below for the simple FITS file interface with CPU Kernel.

  1. A SuperScene that holds global parameter state and parameter bounds, and can be used to check out sub-scenes.

    sceneDB = LinkedSuperScene(sourcecat=cat, bands=bands,
                           statefile=os.path.join(config.outdir, "final_scene.fits"),
                           roi=cat["roi"],
                           bounds_kwargs=bounds_kwargs,
                           target_niter=config.sampling_draws)

  2. A Patcher object that organizes image pixel data and image meta data.

    class Patcher(FITSPatch, CPUPatchMixin):
      pass

patcher = Patcher(fitsfiles=config.image_names,
                  psfstore=config.psfstore,
                  splinedata=config.splinedatafile,
                  sci_ext=1,
                  unc_ext=2,
                  return_residual=True)

  3. Then a loop can start that checks out sub-scenes, finds the relevant pixel data, and constructs an object that can compute posterior probabilities:

    patchID = 0
# Draw scenes until all sources have gotten target number of iterations of HMC
while sceneDB.undone:
    # Draw the sub-scene and associated information
    # A seed of -1 will choose an available scene at random
    region, active, fixed = sceneDB.checkout_region(seed_index=-1)
    bounds, cov = sceneDB.bounds_and_covs(active["source_index"])

    # Collect the pixel data and meta-data
    patcher.build_patch(region, None, allbands=bands)
    # Transfer data to device, subtract fixed sources, set up parameter transforms
    model, q = patcher.prepare_model(active=active, fixed=fixed,
                                    bounds=bounds, shapes=sceneDB.shape_cols)

  4. Within the loop we will either do optimization or HMC sampling, and then check the scene back in. Here we do sampling using littlemcmc:

        # run HMC, with warmup
    out, step, stats = run_lmc(model, q.copy(),
                              n_draws=config.sampling_draws,
                              warmup=config.warmup,
                              z_cov=cov, full=True,
                              weight=max(10, active["n_iter"].min()),
                              discard_tuned_samples=True,
                              max_treedepth=config.max_treedepth,
                              progressbar=config.progressbar)

    # Add additional information to the output
    final, covs = out.fill(region, active, fixed, model, bounds=bounds,
                          step=step, stats=stats, patchID=patchID)
    # Write results to disk
    write_to_disk(out, config.outroot, model, config)
    # Check in the scene with new parameter values
    sceneDB.checkin_region(final, fixed, config.sampling_draws,
                          block_covs=covs, taskID=0)
    # write current global parameter state to disk as a failsafe
    sceneDB.writeout()
    # increment patch number
    patchID += 1

All Steps#

Stepping back a bit, one might want to do an initial round of optimization of the entire catalog, and then use that as initialization for a sampling phase. The steps to do such a full run might look like the following

  1. Create PSF mixtures for mosaic and/or individual exposures

  2. (optional) Pre-Process (preprocess.py)

    This creates the HDF5 storage files for pixel and meta-data. If slopes are present, make separate stores for mosaic and slope pixels.

  3. Make catalog of initial peaks in forcepho format

    The initial peaks to be fit must be supplied in a FITS binary table in the appropriate format. See ./inputs.rst for details of this format.

  4. Background subtraction & optimization loop

    • Optimize sources in the catalog (optimize.py) using mosaic data with a S/N cap.

    • (optional) Fit for a residual background (background.py) in the mosaic. If significant, put resulting tweak values in config file.

    • (optional) Look for objects missing in the initial catalog.

    • (optional) Re-optimize sources in the catalog (optimize.py) based on mosaic.

    • Replace initialization catalog with the optimization results, including flux uncertainty estimates This is done with

      postprocess.py --root output/<run_name> --catname postop_catalog.fits --mode postop

  5. Sample posterior for source properties (sampling.py).

  6. Post-process to create residual images (if available), show chains, etc…

Contents