ProFit: PSF Fitting Example

Prepare the test data

Next we load a table of data describing GAMA galaxies:

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

##    CATAID sersic.xcen1 sersic.ycen1 sersic.mag1 sersic.mag2 sersic.re1
## 1  265769     122.8069      87.3328    17.08907    17.08907    8.71420
## 2  265911      84.8832      94.5951    16.83217    16.83217    7.05740
## 3  265940      80.3975      52.2173    18.04857    18.04857    5.96290
## 4  265943      53.7013      83.6956    17.70547    17.70547    6.42060
## 5  265981     148.8727     106.4069    16.85617    16.85617   11.00100
## 6  265986      78.1780      79.9999    18.97747    18.97747    6.80575
## 7  265985      51.3467      73.0023    19.03897    19.03897    5.77065
## 8  266033      66.0826      79.0872    17.39137    17.39137    6.25480
## 9  266035     105.1495      78.5289    18.18017    18.18017   11.38780
## 10 266105     264.8947     219.5667    15.45027    15.45027   36.68265
##    sersic.re2 sersic.nser1 sersic.ang2 sersic.axrat2
## 1     17.4284       3.4292    105.8498        0.1541
## 2     14.1148       4.3776    140.8191        0.4891
## 3     11.9258       4.6010    112.2746        0.4875
## 4     12.8412       4.5086    168.2385        0.4179
## 5     22.0020       2.9434     59.6690        0.2621
## 6     13.6115       2.0082    103.2134        0.0423
## 7     11.5413       2.7402      0.7028        0.5400
## 8     12.5096       4.9178     39.5336        0.3336
## 9     22.7756       2.9550     56.3237        0.6222
## 10    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.