ProFit: PSF Fitting Example
In ProFit: Fit Projected 2D Profiles to Galaxy Images

Get the latest version of ProFit:

library(devtools)
install_github('ICRAR/ProFit')

First load the libraries we need:

library(ProFit)

Prepare the test data

Next we load a table of data describing GAMA galaxies:

data('ExampleInit', package="ProFit")
head(ExampleInit, 10)

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

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

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

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.