Quickstart guide

imcascade is a method to fit sources in astronomical images. This is accomplished using the multi-Gaussian expansion method which models the galaxy as a mixture of Gaussians. For more details please read the in-depth example. What follows is a (very) brief introduction to the basic usage of imcascade.

In this short example we will fit an analytic, circular, Sersic profile with \(n = 1.5\), \(r_e = 6\) and total flux, \(F = 250\), we have convolved the profile with a Moffat PSF with \(\alpha = 3\) and \(\gamma = 3\) and added purely Gaussian noise.

In a hidden cell I have intialized the cutout in the 2D array sci and the pixelized PSF saved in the 2D array psf. Below I show to 2D images of each

fig,(ax1,ax2) = plt.subplots(1,2, figsize = (12,6))

ax1.imshow(sci)
ax1.set_title('Image to be fit')

ax2.imshow(psf)
ax2.set_title('PSF image')

plt.show()

The initialize_fitter function is designed to take a pixelized science images (as 2D arrays or fits files) and help initalize a Fitter instance which will be used below to fit galaxies. This function is designed to help users get started using our experiences to help guide some of the decisions, which may not be applicable in all situations. For more details about these choices or other options please see the in depth example for a longer discussion about all possibilities.

from imcascade.fitter import initialize_fitter

fitter = initialize_fitter(sci,psf)

2022-07-12 17:48:19,219 - Fit PSF with 4 components
2022-07-12 17:48:19,220 - Widths: 1.2,1.81,5.28,2.93
2022-07-12 17:48:19,221 - Fluxes: 0.27,0.51,0.03,0.19
2022-07-12 17:48:19,231 - Using 9 components with logarithmically spaced widths to fit galaxy
2022-07-12 17:48:19,232 - 0.91, 1.52, 2.55, 4.27, 7.15, 11.97, 20.04, 33.54, 56.15
2022-07-12 17:48:19,233 - No mask was given, derriving one using sep
2022-07-12 17:48:19,246 - Using sep rms map to calculate pixel weights

This function uses the psf_fitter module to fit the given pixelized psf with a sum of Gaussians, which is required for our method. Next it estimates the effective radius and uses nine logarithmically spaced widths for the Gaussian components ranging from the PSF HWHM to \(10\times r_e\). It then derrives pixel weights and masks using sep (or the gain,exposure time and readnoise to calculate the rms). There are also options to use pre-calculated version of these if the user has them.

Now we will run our least-squares minimization routine

opt_param = fitter.run_ls_min()
print (opt_param)

2022-07-12 17:48:21,366 - Running least squares minimization
2022-07-12 17:48:47,347 - Finished least squares minimization

[ 7.60051373e+01  7.50005090e+01  6.87265676e-01  1.57715338e+00
  1.42502826e+00 -2.57490313e+00  1.53411116e+00  1.81786238e+00
  1.87217228e+00  1.63625068e+00  7.22478506e-01 -1.31211431e+00
  4.65480859e-01 -2.74388677e-05 -9.73687762e-05  1.91394558e-04]

We have printed out the parameters the desribe the best fit model. The first four are the structural parameters (x position, y position, axis ratio and position angle) then the next nine represent the fluxes or weights of the nine components (in logarithmic space) and the final three parameters for the tilted-plane sky model.

Obviously this is non-trivial to parse, which is why we will use our results module and the ImcascadeResults class to help us analyze the results

from imcascade.results import ImcascadeResults
res_class = ImcascadeResults(fitter)
res_class.run_basic_analysis()

{'flux': 252.58989931777214,
 'r20': 2.305274076411796,
 'r50': 5.940674533616872,
 'r80': 12.340978912963354,
 'r90': 17.43337150701717,
 'C80_20': 5.353367323755417,
 'C90_50': 2.934577783779577}

The function .run_basic_analysis() calculates some basic morphological quantities like the total flux and half-light radius. We can see that the best fit parameters match the inputs pretty well!

Another very useful function is the .make_diagnostic_fig() this makes a figure which helps inspect the fit

fig = res_class.make_diagnostic_fig()

/mnt/c/Users/timbl/Documents/files/research/packages/imcascade/imcascade/results.py:707: RuntimeWarning: divide by zero encountered in true_divide
  rms_med = np.median(1./np.sqrt(fitter.weight) )

This makes it easier to see if the fit went catastrophically wrong. This fit looks pretty good! Some examples of issues are if the curve-of-growth does not converge or if there are systematic issues in the residuals. To remedy this one could try using a different set of widths for the components, altering the inital guesses or making sure all the input data is correct.

Next we will explore the Posterior distribution using Dynesty. Specifically we will use the ‘express’ method which uses pre-rendered images to help speed up the run time. The code also automatically checks to see if the package jax is installed. This additionally helps to speed up the computation if availible.

post = fitter.run_dynesty(method = 'express')

Now if we re-initialize the results class, we can calculate uncertainties on the morphological values.

res_class_w_post = ImcascadeResults(fitter)
res_class_w_post.run_basic_analysis()

{'flux': array([249.31304659,   1.36354565,   1.27251959]),
 'r20': array([2.3907311 , 0.01458799, 0.01400204]),
 'r50': array([5.88510278, 0.03870749, 0.03532017]),
 'r80': array([12.00401554,  0.13302446,  0.1238212 ]),
 'r90': array([16.70281417,  0.28503741,  0.2769434 ]),
 'C80_20': array([5.0199006 , 0.03832022, 0.03772956]),
 'C90_50': array([2.83774779, 0.03312112, 0.03330771])}