Example: Surface brightness profile extraction and fitting¶
This thread shows how to read data, extract surface brightness profiles, fit data and extract density profiles with PyProffit. The steps presented here can be replicated using the test data available in the validation folder of the PyProffit package.
Reading data¶
We start by loading the packages:
[1]:
import numpy as np
import pyproffit
import matplotlib.pyplot as plt
Then we can move to the data directory with the os package.
[2]:
import os
# Change this to the proper directory containing your run
os.chdir('../../validation/')
os.listdir()
[2]:
['test_sb.fits',
'gsb.fits',
'expose_mask_37.fits.gz',
'test_density.pdf',
'commands.xcm',
'test_outmod.fits',
'test_script.py',
'.ipynb_checkpoints',
'test_plot_fit.pdf',
'test_rec_stan.pdf',
'reference_depr.fits',
'test_save_all.fits',
'pspcb_gain2_256.rsp',
'b_37.fits.gz',
'reference_psf.dat',
'test_dmfilth.fits',
'Untitled1.ipynb',
'ign **-0.42',
'back_37.fits.gz',
'reference_pymc3.dat',
'comp_rec.pdf',
'test_region.reg',
'test_mgas.pdf',
'mybeta_GP.stan',
'reference_OP.dat',
'pspcb_gain2_256.fak',
'Untitled.ipynb',
'sim.txt',
'lumin.txt']
Now we load the data inside a Data object in PyProffit structure:
[3]:
dat=pyproffit.Data(imglink='b_37.fits.gz',explink='expose_mask_37.fits.gz',
bkglink='back_37.fits.gz')
WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / Equatorial system reference
the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]
- Here imglink=’b_37.fits.gz’_ is the link to the image file (count map) to be loaded.
- The option explink=’expose_mask_37.fits.gz’ allows the user to load an exposure map for vignetting correction. In case this option is left blank, a uniform exposure of 1s is assumed for the observation.
- The option bkglink=’back_37.fits.gz’ allows to load an external background map, which will be used when extracting surface brightness profiles.
The images are then loaded into the Data structure and can be easily accessed as below:
[4]:
fig = plt.figure(figsize=(20,20))
s1=plt.subplot(221)
plt.imshow(np.log10(dat.img),aspect='auto')
s2=plt.subplot(222)
plt.imshow(dat.exposure,aspect='auto')
<ipython-input-4-00e9545509bd>:3: RuntimeWarning: divide by zero encountered in log10
plt.imshow(np.log10(dat.img),aspect='auto')
[4]:
<matplotlib.image.AxesImage at 0x7f5cf4821460>
All the areas with zero exposure will be automatically excluded. We can ignore additional regions using the region method of the Data class, which loads a DS9 region file (in image or FK5 format):
[5]:
dat.region('test_region.reg')
Excluded 2 sources
The exposure in the requested areas has been set to 0. Let’s look at the output:
[6]:
plt.clf()
fig = plt.figure(figsize=(10,10))
plt.imshow(dat.exposure,aspect='auto')
[6]:
<matplotlib.image.AxesImage at 0x7f5cf3f98d00>
<Figure size 432x288 with 0 Axes>
The Data structure also contains the dmfilth method, which can be used to fill the masked areas. The method computes a 2D spline interpolation in between the gaps and generates a Poisson realization of the spline interpolated data, such that the filled holes have similar statistical properties to their surroundings
[7]:
dat.dmfilth()
Applying high-pass filter
Interpolating in the masked regions
Filling holes
[8]:
plt.clf()
fig = plt.figure(figsize=(10,10))
plt.imshow(np.log10(dat.filth),aspect='auto')
<ipython-input-8-010fe1205bf9>:3: RuntimeWarning: divide by zero encountered in log10
plt.imshow(np.log10(dat.filth),aspect='auto')
[8]:
<matplotlib.image.AxesImage at 0x7f5cf471eeb0>
<Figure size 432x288 with 0 Axes>
The image produced by dmfilth is to be compared with the raw image shown above; it is apparent that the sources have been removed and their area has been replaced by a Poisson realization of their interpolated surroundings.
In case a dmfilth image has been generated, the computation of the image centroid and/or of the surface brightness peak to compute the center of the surface brightness profile is done on the dmfilth image rather than on the original image.
Profile extraction¶
Now we define a Profile object in the following way:
[9]:
prof=pyproffit.Profile(dat,center_choice='centroid',maxrad=45.,binsize=20.,centroid_region=30.)
Computing centroid and ellipse parameters using principal component analysis
No approximate center provided, will search for the centroid within a radius of 30 arcmin from the center of the image
Denoising image...
Running PCA...
Centroid position: 272.6425582785415 277.47130902570234
Corresponding FK5 coordinates: 55.71693918958471 -53.642873306648674
Ellipse axis ratio and position angle: 1.112280720144629 -144.82405288826155
Profile class options
The class Profile is designed to contain all the Proffit profile extraction features (not all of them have been implemented yet). The “center_choice” argument specifies the choice of the center:
- center_choice=’centroid’: compute image centroid and ellipticity
- center_choice=’peak’: use brightness peak
- center_choice=’custom_fk5’: use custom center in FK5 coordinates (degrees), provided by the “center_ra” and “center_dec” arguments
- center_choice=’custom_ima’: like custom_fk but with input coordinates in image pixels
The other arguments are the following:
- maxrad: define the maximum radius of the profile (in arcmin)
- binsize: the width of the bins (in arcsec)
- center_ra, center_dec: position of the center (if center_choice=’custom_fk5’ or ‘custom_ima’)
- binsize: minimum bin size in arcsec
- binning=: specify binnig scheme: ‘linear’ (default), ‘log’, or ‘custom’. In the ‘custom’ case, an array with the binning definition should be provided through the option bins=array
- centroid_region: for centroid calculation (center_choice=’centroid’), optionally provide a radius within which the centroid will be computed, instead of the entire image.
Now let’s extract the profile…
[10]:
prof.SBprofile(ellipse_ratio=prof.ellratio,rotation_angle=prof.ellangle+180.)
Here we have extracted a profile in elliptical annuli centered on the image centroid (see above), with an ellipse axis ratio (major/minor) and position angle calculated with principal component analysis. If ellipse_ratio and ellipse_angle are left blank circular annuli are used.
- ellipse_ratio: the ratio of major to minor axis (a/b) of the ellipse (default=1, i.e. circular annuli)
- rotation_angle: rotation angle of the ellipse from the R.A. axis (default=0)
- angle_low, angle_high: in case of profile extraction in sectors, the position angle of the minimum and maximum angles of the sector, with 0 equivalent to the R.A. axis (default=None, i.e. the entire azimuth)
Now let’s plot the profile…
[11]:
prof.Plot()
<Figure size 432x288 with 0 Axes>
Defining a model¶
Models can be defined using the Model class. PyProffit includes several popular built-in models, however the Model structure is designed to be compatible with any custom Python function (see below)
[12]:
mod=pyproffit.Model(pyproffit.BetaModel)
To check the parameters of the BetaModel function,
[13]:
mod.parnames
[13]:
('beta', 'rc', 'norm', 'bkg')
Any user-defined Python function operating on NumPy arrays can be defined here, see below.
Fitting the data¶
To fit the extracted profiles PyProffit provides the Fitter class, which takes a Profile and a Model object as input:
[14]:
fitobj=pyproffit.Fitter(model=mod, profile=prof, beta=0.7, rc=2.,norm=-2,bkg=-4)
fitobj.Migrad()
┌──────────────────────────────────┬──────────────────────────────────────┐
│ FCN = 143.5 │ Nfcn = 273 │
│ EDM = 1.13e-05 (Goal: 0.0002) │ │
├───────────────┬──────────────────┼──────────────────────────────────────┤
│ Valid Minimum │ Valid Parameters │ No Parameters at limit │
├───────────────┴──────────────────┼──────────────────────────────────────┤
│ Below EDM threshold (goal x 10) │ Below call limit │
├───────────────┬──────────────────┼───────────┬─────────────┬────────────┤
│ Covariance │ Hesse ok │ Accurate │ Pos. def. │ Not forced │
└───────────────┴──────────────────┴───────────┴─────────────┴────────────┘
┌───┬──────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐
│ │ Name │ Value │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit- │ Limit+ │ Fixed │
├───┼──────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤
│ 0 │ beta │ 0.672 │ 0.016 │ │ │ │ │ │
│ 1 │ rc │ 3.27 │ 0.15 │ │ │ │ │ │
│ 2 │ norm │ -1.412 │ 0.016 │ │ │ │ │ │
│ 3 │ bkg │ -3.733 │ 0.020 │ │ │ │ │ │
└───┴──────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘
┌──────┬─────────────────────────────────────────┐
│ │ beta rc norm bkg │
├──────┼─────────────────────────────────────────┤
│ beta │ 0.000267 0.00226 -0.000138 0.000215 │
│ rc │ 0.00226 0.0222 -0.00184 0.00158 │
│ norm │ -0.000138 -0.00184 0.000258 -7.96e-05 │
│ bkg │ 0.000215 0.00158 -7.96e-05 0.000386 │
└──────┴─────────────────────────────────────────┘
Best fit chi-squared: 143.525 for 135 bins and 131 d.o.f
Reduced chi-squared: 1.09561
The fit provides the best-fit parameters with their errors and the covariance matrix. The fit statistic is provided at the bottom. All the results of the optimization can be accessed through the fitobj.out attribute of the Fitter class, which contains the output of the iminuit migrad command.
Now we can plot the data together with the best fitting model
[15]:
prof.Plot(model=mod)
<Figure size 432x288 with 0 Axes>
That’s nice; now if instead of :math:chi^2` <https://pyproffit.readthedocs.io/en/latest/pyproffit.html#pyproffit.fitting.ChiSquared>`__ we want to fit the counts with C statistic
[16]:
fitobj=pyproffit.Fitter(model=mod, method='cstat', profile=prof, beta=0.7, rc=2.,norm=-2, bkg=-4, fitlow=0., fithigh=30.)
fitobj.Migrad()
┌──────────────────────────────────┬──────────────────────────────────────┐
│ FCN = 86.67 │ Nfcn = 316 │
│ EDM = 9.35e-06 (Goal: 0.0002) │ │
├───────────────┬──────────────────┼──────────────────────────────────────┤
│ Valid Minimum │ Valid Parameters │ No Parameters at limit │
├───────────────┴──────────────────┼──────────────────────────────────────┤
│ Below EDM threshold (goal x 10) │ Below call limit │
├───────────────┬──────────────────┼───────────┬─────────────┬────────────┤
│ Covariance │ Hesse ok │ Accurate │ Pos. def. │ Not forced │
└───────────────┴──────────────────┴───────────┴─────────────┴────────────┘
┌───┬──────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐
│ │ Name │ Value │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit- │ Limit+ │ Fixed │
├───┼──────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤
│ 0 │ beta │ 0.658 │ 0.019 │ │ │ │ │ │
│ 1 │ rc │ 3.17 │ 0.17 │ │ │ │ │ │
│ 2 │ norm │ -1.406 │ 0.017 │ │ │ │ │ │
│ 3 │ bkg │ -3.75 │ 0.04 │ │ │ │ │ │
└───┴──────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘
┌──────┬─────────────────────────────────────────┐
│ │ beta rc norm bkg │
├──────┼─────────────────────────────────────────┤
│ beta │ 0.000372 0.00298 -0.000178 0.000673 │
│ rc │ 0.00298 0.0273 -0.00216 0.00474 │
│ norm │ -0.000178 -0.00216 0.000284 -0.000235 │
│ bkg │ 0.000673 0.00474 -0.000235 0.00191 │
└──────┴─────────────────────────────────────────┘
Best fit C-statistic: 86.6662 for 135 bins and 131 d.o.f
Reduced C-statistic: 0.661574
Here we have restricted the fitting range to be between 0 and 30 arcmin through the fitlow and fithigh arguments. Note that the “reduced C-statistic” cannot be used directly as a goodness-of-fit indicator.
[17]:
prof.Plot(model=mod)
<Figure size 432x288 with 0 Axes>
Inspecting the results¶
After running Migrad the Fitter object contains a minuit object which can be used to run all iminuit tools. For instance, we can run minos to get more accurate errors on the parameters:
[18]:
fitobj.minuit.minos()
[18]:
| FCN = 86.67 | Nfcn = 673 | |||
| EDM = 9.35e-06 (Goal: 0.0002) | ||||
| Valid Minimum | Valid Parameters | No Parameters at limit | ||
| Below EDM threshold (goal x 10) | Below call limit | |||
| Covariance | Hesse ok | Accurate | Pos. def. | Not forced |
| Name | Value | Hesse Error | Minos Error- | Minos Error+ | Limit- | Limit+ | Fixed | |
|---|---|---|---|---|---|---|---|---|
| 0 | beta | 0.658 | 0.019 | -0.019 | 0.020 | |||
| 1 | rc | 3.17 | 0.17 | -0.16 | 0.17 | |||
| 2 | norm | -1.406 | 0.017 | -0.017 | 0.017 | |||
| 3 | bkg | -3.75 | 0.04 | -0.05 | 0.04 |
| beta | rc | norm | bkg | |||||
|---|---|---|---|---|---|---|---|---|
| Error | -0.019 | 0.020 | -0.16 | 0.17 | -0.017 | 0.017 | -0.05 | 0.04 |
| Valid | True | True | True | True | True | True | True | True |
| At Limit | False | False | False | False | False | False | False | False |
| Max FCN | False | False | False | False | False | False | False | False |
| New Min | False | False | False | False | False | False | False | False |
| beta | rc | norm | bkg | |
|---|---|---|---|---|
| beta | 0.000372 | 0.00298 (0.936) | -0.000178 (-0.548) | 0.000673 (0.797) |
| rc | 0.00298 (0.936) | 0.0273 | -0.00216 (-0.774) | 0.00474 (0.656) |
| norm | -0.000178 (-0.548) | -0.00216 (-0.774) | 0.000284 | -0.000235 (-0.319) |
| bkg | 0.000673 (0.797) | 0.00474 (0.656) | -0.000235 (-0.319) | 0.00191 |
The draw_mncontour method allows to compute parameter covariance plots
[19]:
fitobj.minuit.draw_mncontour('beta', 'rc', cl=(0.68, 0.9, 0.99))
[19]:
<matplotlib.contour.ContourSet at 0x7f5cf418c0a0>
Likelihood profiles for each parameter can be created using the draw_mnprofile method,
[20]:
rcval, fcnval = fitobj.minuit.draw_mnprofile('rc')
Setting constraints on the parameters¶
To help convergence or avoid bad results it is often useful to restrict the range of acceptable values for a parameter or fix its value. This can be done easily by accessing the minuit attribute of the Fitter class.
For instance, let’s assume we want to fix the value of the beta parameter to 0.7 and refit
[21]:
fitobj.minuit.fixed['beta'] = True
fitobj.minuit.migrad()
[21]:
| FCN = 86.67 | Nfcn = 15613 | |||
| EDM = 1.38e-13 (Goal: 0.0002) | ||||
| Valid Minimum | Valid Parameters | No Parameters at limit | ||
| Below EDM threshold (goal x 10) | Below call limit | |||
| Covariance | Hesse ok | Accurate | Pos. def. | Not forced |
| Name | Value | Hesse Error | Minos Error- | Minos Error+ | Limit- | Limit+ | Fixed | |
|---|---|---|---|---|---|---|---|---|
| 0 | beta | 0.658 | 0.019 | -0.019 | 0.020 | yes | ||
| 1 | rc | 3.17 | 0.06 | -0.16 | 0.17 | |||
| 2 | norm | -1.406 | 0.014 | -0.017 | 0.017 | |||
| 3 | bkg | -3.750 | 0.026 | -0.047 | 0.041 |
| beta | rc | norm | bkg | |||||
|---|---|---|---|---|---|---|---|---|
| Error | -0.019 | 0.020 | -0.16 | 0.17 | -0.017 | 0.017 | -0.05 | 0.04 |
| Valid | True | True | True | True | True | True | True | True |
| At Limit | False | False | False | False | False | False | False | False |
| Max FCN | False | False | False | False | False | False | False | False |
| New Min | False | False | False | False | False | False | False | False |
| beta | rc | norm | bkg | |
|---|---|---|---|---|
| beta | 0 | 0 | 0 | 0 |
| rc | 0 | 0.00338 | -0.000729 (-0.889) | -0.000651 (-0.424) |
| norm | 0 | -0.000729 (-0.889) | 0.000199 | 8.67e-05 (0.233) |
| bkg | 0 | -0.000651 (-0.424) | 8.67e-05 (0.233) | 0.000698 |
If on top of that we want to restrict parameter rc to lie within the range (0,10):
[22]:
fitobj.minuit.limits['rc'] = (0, 10)
fitobj.minuit.migrad()
[22]:
| FCN = 86.67 | Nfcn = 15649 | |||
| EDM = 6.24e-13 (Goal: 0.0002) | ||||
| Valid Minimum | Valid Parameters | No Parameters at limit | ||
| Below EDM threshold (goal x 10) | Below call limit | |||
| Covariance | Hesse ok | Accurate | Pos. def. | Not forced |
| Name | Value | Hesse Error | Minos Error- | Minos Error+ | Limit- | Limit+ | Fixed | |
|---|---|---|---|---|---|---|---|---|
| 0 | beta | 0.658 | 0.019 | -0.019 | 0.020 | yes | ||
| 1 | rc | 3.17 | 0.06 | -0.16 | 0.17 | 0 | 10 | |
| 2 | norm | -1.406 | 0.014 | -0.017 | 0.017 | |||
| 3 | bkg | -3.750 | 0.026 | -0.047 | 0.041 |
| beta | rc | norm | bkg | |||||
|---|---|---|---|---|---|---|---|---|
| Error | -0.019 | 0.020 | -0.16 | 0.17 | -0.017 | 0.017 | -0.05 | 0.04 |
| Valid | True | True | True | True | True | True | True | True |
| At Limit | False | False | False | False | False | False | False | False |
| Max FCN | False | False | False | False | False | False | False | False |
| New Min | False | False | False | False | False | False | False | False |
| beta | rc | norm | bkg | |
|---|---|---|---|---|
| beta | 0 | 0 | 0 | 0 |
| rc | 0 | 0.00338 | -0.000729 (-0.889) | -0.000651 (-0.424) |
| norm | 0 | -0.000729 (-0.889) | 0.000199 | 8.67e-05 (0.233) |
| bkg | 0 | -0.000651 (-0.424) | 8.67e-05 (0.233) | 0.000698 |
Subtracting the background¶
Now if we want to re-extract the profile with a different binning and subtract the background, we can create a new Profile object…
[23]:
p2=pyproffit.Profile(dat,center_choice='custom_ima',center_ra=prof.cx,center_dec=prof.cy,
maxrad=25.,binsize=20.,binning='log')
p2.SBprofile(ellipse_ratio=prof.ellratio, rotation_angle=prof.ellangle)
Corresponding FK5 coordinates: 55.72393261623938 -53.647031714234124
[24]:
p2.Plot()
<Figure size 432x288 with 0 Axes>
We will use the results of the previous fit stored in the Fitter object to subtract the background. The Backsub method of the Profile class reads the bkg parameter, subtracts its value from the brightness profile and adds the statistical error in quadrature:
[25]:
p2.Backsub(fitobj)
p2.Plot()
<Figure size 432x288 with 0 Axes>
Fitting with a custom function¶
To fit the data with any custom function we can simply define a Python function with the desired model and pass it to a PyProffit Model object, with the following structure:
[26]:
# Create your own model, here some sort of cuspy beta model
def myModel(x,beta,rc,alpha,norm,bkg):
term1 = np.power(x/rc,-alpha) * np.power(1. + (x/rc) ** 2, -3 * beta + alpha/2.)
n2 = np.power(10., norm)
b2 = np.power(10., bkg)
return n2 * term1 + b2
[27]:
# Pass the user-defined function to a new Model structure
cuspmod = pyproffit.Model(myModel)
Now we fit the model to the data, as done previously.
We can also fix the value of some of the parameters to help the convergence. This is done easily by modifying the iminuit.fixed variable for the parameter, as shown above.
[28]:
fitcusp = pyproffit.Fitter(model=cuspmod, profile=prof, beta=0.35, rc=1., alpha=0., norm=-1.7, bkg=-3.8)
fitcusp.Migrad()
┌──────────────────────────────────┬──────────────────────────────────────┐
│ FCN = 140.6 │ Nfcn = 407 │
│ EDM = 3.49e-05 (Goal: 0.0002) │ │
├───────────────┬──────────────────┼──────────────────────────────────────┤
│ Valid Minimum │ Valid Parameters │ No Parameters at limit │
├───────────────┴──────────────────┼──────────────────────────────────────┤
│ Below EDM threshold (goal x 10) │ Below call limit │
├───────────────┬──────────────────┼───────────┬─────────────┬────────────┤
│ Covariance │ Hesse ok │ Accurate │ Pos. def. │ Not forced │
└───────────────┴──────────────────┴───────────┴─────────────┴────────────┘
┌───┬───────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐
│ │ Name │ Value │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit- │ Limit+ │ Fixed │
├───┼───────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤
│ 0 │ beta │ 0.490 │ 0.017 │ │ │ │ │ │
│ 1 │ rc │ 2.95 │ 0.23 │ │ │ │ │ │
│ 2 │ alpha │ -0.14 │ 0.09 │ │ │ │ │ │
│ 3 │ norm │ -1.33 │ 0.05 │ │ │ │ │ │
│ 4 │ bkg │ -3.742 │ 0.021 │ │ │ │ │ │
└───┴───────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘
┌───────┬───────────────────────────────────────────────────┐
│ │ beta rc alpha norm bkg │
├───────┼───────────────────────────────────────────────────┤
│ beta │ 0.000303 0.00353 0.000822 -0.00054 0.000253 │
│ rc │ 0.00353 0.0533 0.0172 -0.0105 0.00247 │
│ alpha │ 0.000822 0.0172 0.00819 -0.0043 0.000495 │
│ norm │ -0.00054 -0.0105 -0.0043 0.0025 -0.000328 │
│ bkg │ 0.000253 0.00247 0.000495 -0.000328 0.000438 │
└───────┴───────────────────────────────────────────────────┘
Best fit chi-squared: 140.643 for 135 bins and 130 d.o.f
Reduced chi-squared: 1.08187
[29]:
prof.Plot(model=cuspmod)
<Figure size 432x288 with 0 Axes>
