RDI with a PSF Library

pyKLIP supports RDI with a PSF Library using the pyklip.rdi.PSFLibrary class. The PSF Library class holds a correlation matrix of all frames with each other, so it knows which frames are good reference frames for a certain frame. This means that all the data, science data and reference data, should be in the PSF libary. In this sense, it is good to generate the PSF Library class once, spend the time to compute the correlation matrix once, and then save the correlation to file using the save_correlation() function inside of PSFLibrary. This way, the PSF library can be applied on multiple science targets, and the correlation matrix needs to just be computed once.

Set up a PSF Library

Here we will assume you have a 3-D cube of frames from a long series of data (e.g. a night of observing, a full survey) that is already in the variable psflib_imgs. This could be made from your own code, or from dataset.input after you read in a large list of files into the Data object. All of these images need to be aligned to some defined aligned_center, a [x,y] array. If the images haven’t been aligned to this common center already, you can use pyklip.klip.align_and_scale() to register each frame. If you are using GPI IFS data, you can also use the pyklip.instruments.Instrument.Data.spectral_collapse() to simultaeously center and collapse a dataset. You will also need psflib_filenames, an array of filenames, so that each frame has a corresponding filename. This again can be taken from dataset.filenames if you don’t have it already, but have read in all the files into a Data object. You want to make sure that the filenames are accurate, as this is how we figure out which frames to exclude for RDI (i.e. we don’t want to use data in the science sequence as reference images in a RDI reduction).

Now here we demonstrate how to make the PSFLibrary object, compute the correlation matrix, and save the correlation matrix. Note that with a lot of files, generating the correlation matrix can take a long time.

import pyklip.rdi as rdi

# make the PSF library
# we need to compute the correlation matrix of all images vs each other since we haven't computed it before
psflib = rdi.PSFLibrary(psflib_imgs, aligned_center, psflib_filenames, compute_correlation=True)

# save the correlation matrix to disk so that we also don't need to recomptue this ever again
# In the future we can just pass in the correlation matrix into the PSFLibrary object rather than having it compute it
psflib.save_correlation("corr_matrix.fits", overwrite=True)

Then, in the future, we don’t need to recompute the correlation matrix when setting up the PSF library. Instead we can read it in and regenerate the PSFLibrary quickly.

import pyklip.rdi as rdi
import astropy.io.fits as fits

# read in the correlation matrix we already saved
corr_matrix_hdulist = fits.open("corr_matrix.fits")
corr_matrix = corr_matrix_hdulist[0].data

# make the PSF library again, this time we have the correlation matrix
psflib = rdi.PSFLibrary(psflib_imgs, aligned_center, psflib_filenames, correlation_matrix=corr_matrix)

Running RDI on a dataset

Now let’s assume you have a dataset, a object that implements pyklip.instruments.Instrument.Data. The input files of dataset are also part of the PSF library (e.g. part of the data taken in that night, or as part of the survey). We will then prepare the PSF library for this dataset. This basically invalidates the frames in the PSF library that belong to this target so they aren’t used as reference images, so there won’t be any self-subtraction. This is where we use the filenames to match frames, so it is important that dataset.filenames matches the filenames passed into the PSF library.

Then, it’s as simple as running KLIP. We pass in the same aligned_center that all the images in the PSF library is aligned to so that this reduction also aligns the science images to that center.

# now we need to prepare the PSF library to reduce this dataset
# what this does is tell the PSF library to not use files from this star in the PSF library for RDI
psflib.prepare_library(dataset)

# now we can run RDI klip
# as all RDI images are aligned to aligned_center, we need to pass in that aligned_center into KLIP
numbasis=[1,5,10,20,50] # number of KL basis vectors to use to model the PSF. We will try several different ones
maxnumbasis=150 # maximum number of most correlated PSFs to do PCA reconstruction with
annuli=3
subsections=4 # break each annulus into 4 sectors
parallelized.klip_dataset(dataset, outputdir="data/", fileprefix="pyklip_k150a3s4m1", annuli=annuli,
                        subsections=subsections, numbasis=numbasis, maxnumbasis=maxnumbasis, mode="RDI",
                        aligned_center=aligned_center, psf_library=psflib, movement=1)

RDI Forward modeling

You can also use RDI with most of the forward modeling features in pyKLIP. The following forward modeling features can be used with RDI with slight modifications:

Foward-Model Astrometry and Photometry

Forward Model Matched Filter (FMMF)

Disk Foward Modelling (DiskFM)

The main changes to make are just to set mode='RDI' and include the keyword psf_library=psflib in pyklip.fm.klip_dataset(). The rest will be taken care of under the hood. It is also possible to combine RDI with ADI/SDI.

Note that the current PSF library does not distinguish spectral channels so it is not well suited for spectral RDI capabilities, and thus, the forward model spectral extraction is not currently compatible with RDI.