ProFit: PSF Fitting Example

Aaron Robotham

2017-12-08

Prepare the test data

First load the libraries we need:

library(knitr)
library(ProFit)
## Loading required package: Rcpp
## Warning: package 'Rcpp' was built under R version 3.4.2
## Loading required package: FITSio
library(FITSio)

Next we load a table of data describing GAMA galaxies:

data('ExampleInit', package="ProFit")
kable(head(ExampleInit, 10))
CATAID sersic.xcen1 sersic.ycen1 sersic.mag1 sersic.mag2 sersic.re1 sersic.re2 sersic.nser1 sersic.ang2 sersic.axrat2
265769 122.8069 87.3328 17.08908 17.08908 8.71420 17.4284 3.4292 105.8498 0.1541
265911 84.8832 94.5951 16.83217 16.83217 7.05740 14.1148 4.3776 140.8191 0.4891
265940 80.3975 52.2173 18.04857 18.04857 5.96290 11.9258 4.6010 112.2746 0.4875
265943 53.7013 83.6956 17.70547 17.70547 6.42060 12.8412 4.5086 168.2385 0.4179
265981 148.8727 106.4069 16.85618 16.85618 11.00100 22.0020 2.9434 59.6690 0.2621
265986 78.1780 79.9999 18.97747 18.97747 6.80575 13.6115 2.0082 103.2134 0.0423
265985 51.3467 73.0023 19.03898 19.03898 5.77065 11.5413 2.7402 0.7028 0.5400
266033 66.0826 79.0872 17.39138 17.39138 6.25480 12.5096 4.9178 39.5336 0.3336
266035 105.1495 78.5289 18.18017 18.18017 11.38780 22.7756 2.9550 56.3237 0.6222
266105 264.8947 219.5667 15.45027 15.45027 36.68265 73.3653 3.9704 94.4409 0.3477

Here we will use an SDSS example:

datasource='SDSS'

Now we can extract out the example files we have available for fitting by checking the contents of the directory containing the example FITS files:

ExampleFiles=list.files(system.file("extdata",datasource,package="ProFit"))
ExampleIDs=unlist(strsplit(ExampleFiles[grep('fitim',ExampleFiles)],'fitim.fits'))
ExampleIDs
##  [1] "G265911" "G265940" "G266033" "G266035" "G266662" "G267199" "G267489"
##  [8] "G267525" "G278109" "G279148"

There are 10 example galaxies included. Here we run example 1:

useID=ExampleIDs[1]
psf=readFITS(system.file("extdata", paste(datasource,'/',useID,'psfim.fits',sep=''),package="ProFit"))$imDat
psf=psf/sum(psf)
psfsigma=sqrt(abs(psf))/200 #To get reasonable PSF errors

We can check the image of the PSF wiht magimage:

temp=magimage(psf,lo=0,hi=1)
contour(temp, add=T, col='red', drawlabels=FALSE)

To check the profile a few 1D plots can be useful:

magplot(psf[13,],type='l')
lines(psf[12,],lty=2,col='red')
lines(psf[14,],lty=2,col='blue')
lines(psf[11,],lty=3,col='red')
lines(psf[15,],lty=3,col='blue')

The red and blue lines fall pretty much on top of each other, so there is not much ellipticity.

Setup the fitting data structures

We can use ProFit to fit an analytic Moffat function to the PSF to properly characterise this.

modellist=list(
  moffat=list(
    xcen=dim(psf)[1]/2,
    ycen=dim(psf)[2]/2,
    mag=0,
    fwhm=2.5,
    con=3,
    ang=0,
    axrat=0.9,
    box=0
  )
)
modellist
## $moffat
## $moffat$xcen
## [1] 12.5
## 
## $moffat$ycen
## [1] 12.5
## 
## $moffat$mag
## [1] 0
## 
## $moffat$fwhm
## [1] 2.5
## 
## $moffat$con
## [1] 3
## 
## $moffat$ang
## [1] 0
## 
## $moffat$axrat
## [1] 0.9
## 
## $moffat$box
## [1] 0

We can check what this default model looks like:

psfmodel=profitMakeModel(modellist,dim=c(25,25))$z
magimage(psfmodel)
magimage(abs(psfmodel-psf)/psfmodel,zlim=c(0,1))

We will fit everything:

tofit=list(
  moffat=list(
    xcen=TRUE,
    ycen=TRUE,
    mag=TRUE,
    fwhm=TRUE,
    con=TRUE,
    ang=TRUE,
    axrat=TRUE,
    box=TRUE
  )
)

And choose sensible options for which parameters to fit in log-space:

tolog=list(
  moffat=list(
    xcen=FALSE,
    ycen=FALSE,
    mag=FALSE,
    fwhm=FALSE,
    con=TRUE,
    ang=FALSE,
    axrat=TRUE,
    box=FALSE
  )
)

The intervals will be generous:

intervals=list(
  moffat=list(
    xcen=list(lim=c(0,25)),
    ycen=list(lim=c(0,25)),
    mag=list(lim=c(-1,1)),
    fwhm=list(lim=c(0.1,10)),
    con=list(lim=c(1,20)),
    ang=list(lim=c(-180,360)),
    axrat=list(lim=c(0.1,1)),
    box=list(lim=c(-1,1))
  )
)

Now setup the Data structure we need for fitting, where we will use Normal likelihoods:

Data=profitSetupData(image=psf, mask=(psf==0), sigma=psfsigma, modellist=modellist, tofit=tofit, tolog=tolog, intervals=intervals, magzero=0, algo.func='optim', verbose=TRUE, like.func='norm')

Check our rough model:

profitLikeModel(parm=Data$init, Data=Data, makeplots=TRUE, plotchisq=TRUE)

The initial rough model for the SDSS does not look great- clearly the PSF is much rounder than we guessed.

To stop the guess-work we can now optimise the model with BFGS:

optimfit=optim(Data$init, profitLikeModel, method='BFGS', Data=Data, control=list(fnscale=-1))

Check the final result:

profitLikeModel(optimfit$par, Data,makeplots=TRUE, plotchisq=TRUE)

You can go even further with a full MCMC fit (I can assure you that in this case it is not worth the effort, but run if you must!):

Data$algo.func="LD"

LDfit=LaplacesDemon(profitLikeModel, Initial.Values=optimfit$par, Data=Data,
  Iterations=1e4, Algorithm='CHARM', Thinning=1, Specs=list(alpha.star=0.44))

Check the final result:

profitLikeModel(LDfit$Summary2[,1], Data,makeplots=TRUE, plotchisq=TRUE)

Since the full MCMC achieves the same fit (more or less) we will procede with the optim fit. In the resultant fit we see that the FWHM is given as ~3, which given the SDSS pixel scale (0.339 asec/pix) is ~1 asec, which is pretty good for SDSS imaging. The PSF is preferred as being close to an axial ratio ~1 (or ~0 in log-space). We do find significant boxiness, so the PSF is not actually perfectly circular. We can see this visually in fact:

PSFmodellist=profitRemakeModellist(optimfit$par, modellist, tofit=tofit, tolog=tolog, intervals=intervals)$modellist
PSFmodellist$moffat$xcen=12.5
PSFmodellist$moffat$ycen=12.5
PSFmodellist$moffat$mag=0
psfmodel=profitMakeModel(PSFmodellist, dim=c(25,25))$z
contour(magimage(psfmodel), add=TRUE, col='red', drawlabels=FALSE)

Notice in the above we set a few parameters to be exactly where we want, e.g. the PSF should be in the middle of our 25x25 image matrix (the fit was at 12.51 and 12.51 rather than 12.5 and 12.5) and the magnitude should be exactly 0 (i.e. the integral of the PSF sums to 1). We now have an analytic means of describing the SDSS PSF. We can use this for subsequent fitting.

Use the analytic PSF to fit a target galaxy

useID=ExampleIDs[1]
image=readFITS(system.file("extdata", paste(datasource,'/',useID,'fitim.fits',sep=''), package="ProFit"))$imDat
sigma=readFITS(system.file("extdata", paste(datasource,'/',useID,'sigma.fits',sep=''), package="ProFit"))$imDat
segim=readFITS(system.file("extdata", paste(datasource,'/',useID,'segim.fits',sep=''), package="ProFit"))$imDat

Next we extract parameters for a very rough model (not meant to look too good yet):

useIDnum=as.integer(strsplit(useID,'G')[[1]][2])
useloc=which(ExampleInit$CATAID==useIDnum)

Setup the fitting data structures

For our initial model we treat component 1 as the putative bulge and component 2 as the putative disk. We are going to attempt a fit where the disk is forced to have nser=1 and the bulge has an axial ratio of 1.

modellist=list(
  sersic=list(
    xcen= c(dim(image)[1]/2, dim(image)[1]/2),
    ycen= c(dim(image)[2]/2, dim(image)[2]/2),
    mag= c(ExampleInit$sersic.mag1[useloc], ExampleInit$sersic.mag2[useloc]),
    re= c(ExampleInit$sersic.re1[useloc], ExampleInit$sersic.re2[useloc])*
      if(datasource=='KiDS'){1}else{0.2/0.339},
    nser= c(ExampleInit$sersic.nser1[useloc], 1),  #Disk is initially nser=1
    ang= c(ExampleInit$sersic.ang2[useloc], ExampleInit$sersic.ang2[useloc]),
    axrat= c(1, ExampleInit$sersic.axrat2[useloc]),  #Bulge is initially axrat=1
    box=c(0, 0)
  )
)
modellist
## $sersic
## $sersic$xcen
## [1] 50.5 50.5
## 
## $sersic$ycen
## [1] 50.5 50.5
## 
## $sersic$mag
## [1] 16.83217 16.83217
## 
## $sersic$re
## [1] 4.163658 8.327316
## 
## $sersic$nser
## [1] 4.3776 1.0000
## 
## $sersic$ang
## [1] 140.8191 140.8191
## 
## $sersic$axrat
## [1] 1.0000 0.4891
## 
## $sersic$box
## [1] 0 0

The pure model (no PSF):

magimage(profitMakeModel(modellist,dim=dim(image)))

The original image:

magimage(image)

The convolved model (with PSF):

magimage(profitMakeModel(modellist,dim=dim(image),psf=psfmodel))

Next we define our list of what we want to fit (where TRUE means we will fit it later):

tofit=list(
  sersic=list(
    xcen= c(TRUE,NA), #We fit for xcen and tie the two togther
    ycen= c(TRUE,NA), #We fit for ycen and tie the two togther
    mag= c(TRUE,TRUE), #Fit for both
    re= c(TRUE,TRUE), #Fit for both
    nser= c(TRUE,FALSE), #Fit for bulge
    ang= c(FALSE,TRUE), #Fit for disk
    axrat= c(FALSE,TRUE), #Fit for disk
    box= c(FALSE,FALSE) #Fit for neither
  )
)

Now we define what parameters should be fitted in log space:

tolog=list(
  sersic=list(
    xcen= c(FALSE,FALSE),
    ycen= c(FALSE,FALSE),
    mag= c(FALSE,FALSE),
    re= c(TRUE,TRUE), #re is best fit in log space
    nser= c(TRUE,TRUE), #nser is best fit in log space
    ang= c(FALSE,FALSE),
    axrat= c(TRUE,TRUE), #axrat is best fit in log space
    box= c(FALSE,FALSE)
  )
)

The hard intervals should also be specified in linear space:

intervals=list(
  sersic=list(
    xcen=list(lim=c(0,300),lim=c(0,300)),
    ycen=list(lim=c(0,300),lim=c(0,300)),
    mag=list(lim=c(10,30),lim=c(10,30)),
    re=list(lim=c(1,100),lim=c(1,100)),
    nser=list(lim=c(0.5,20),lim=c(0.5,20)),
    ang=list(lim=c(-180,360),lim=c(-180,360)),
    axrat=list(lim=c(0.1,1),lim=c(0.1,1)),
    box=list(lim=c(-1,1),lim=c(-1,1))
  )
)

Setup the data structure we need for optimisation, taking a few seconds to find the optimal convolution method:

Data=profitSetupData(image=image, sigma=sigma, segim=segim, psf=psfmodel,
                     modellist=modellist, tofit=tofit, tolog=tolog, intervals=intervals,
                     magzero=0, algo.func='optim', like.func="t", verbose=TRUE)

Do some fitting with the new model PSF:

We will try optim BFGS:

optimfitMod=optim(Data$init, profitLikeModel, method='BFGS', Data=Data,
                  control=list(fnscale=-1))

The best optim BFGS fit is given by:

optimfitMod$par

Check it out:

profitLikeModel(optimfitMod$par,Data,makeplots=TRUE,whichcomponents=list(sersic=1))
profitLikeModel(optimfitMod$par,Data,makeplots=TRUE,whichcomponents=list(sersic=2))
profitLikeModel(optimfitMod$par,Data,makeplots=TRUE,whichcomponents=list(sersic='all'))
modeloptim=profitRemakeModellist(optimfitMod$par,Data$modellist,Data$tofit,Data$tolog)$modellist
profitEllipsePlot(Data,modeloptim,pixscale=0.339,FWHM=1,SBlim=26)

Now we can try the empirical PSF instead (for comparison):

Data=profitSetupData(image=image, sigma=sigma, segim=segim, psf=psf,
  modellist=modellist, tofit=tofit, tolog=tolog, intervals=intervals, magzero=0,
  algo.func='optim', like.func="t", verbose=TRUE)

Do some fitting with the empirical model PSF:

We will try optim BFGS:

optimfitEmp=optim(Data$init, profitLikeModel, method='BFGS', Data=Data,
                  control=list(fnscale=-1))

The best optim L-BFGS-B fit is given by:

optimfitEmp$par

Check it out:

profitLikeModel(optimfitEmp$par,Data,makeplots=TRUE,whichcomponents=list(sersic=1))
profitLikeModel(optimfitEmp$par,Data,makeplots=TRUE,whichcomponents=list(sersic=2))
profitLikeModel(optimfitEmp$par,Data,makeplots=TRUE,whichcomponents=list(sersic='all'))
modeloptim=profitRemakeModellist(optimfitEmp$par,Data$modellist,Data$tofit,Data$tolog)$modellist
profitEllipsePlot(Data,modeloptim,pixscale=0.339,FWHM=1,SBlim=26)

Fitting using the empirical PSF gives similar results and best-fit LL (slightly better LL), but a larger bulge Re and a brighter bulge magnitude. It makes sense the model versus empirical PSF would disagree most for bulge parameters- these are dominated by the core region where the resolution is barely above the PSF itself.