ProFound/ProFit: A Complex Fit
In ProFound: Photometry Tools

Get the latest version of ProFound and ProFit:

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

Set global evaluate (basically TRUE for GitHub version and FALSE for CRAN):

evalglobal=FALSE

First load the libraries we need:

library(ProFound)
library(ProFit)

image=readFITS(system.file("extdata", 'VIKING/mystery_VIKING_Z.fits', package="ProFound"))
profound=profoundProFound(image, magzero=30, verbose=TRUE, plot=TRUE, boundstats=TRUE, rotstats=TRUE)

We can identify the large group of central potentially confused sources from the group output from ProFound:

profound$group$groupsegID[1,]

The main group has 8 segments, which we can view easily:

magimage(profound$group$groupim==1)

We can check the properties of the group 1 segments too:

group1segstats=profound$segstats[profound$segstats$segID %in% unlist(profound$group$groupsegID[1,"segID"]),]
group1segstats

In our Full Monty vignette example we detail how to extract a reasonable PSF. We use the solution here.

psf_modellist=list(
  moffat=list(
    xcen=75/2,
    ycen=75/2,
    mag=0,
    fwhm=3.8,
    con=2.04,
    ang=0,
    axrat=1,
    box=0
  )
)
psf_model=profitMakeModel(modellist=psf_modellist, dim=c(75,75))$z

We can see where this FWHM lies in comparison to the values estimated from ProFound.

maghist(profound$segstats$R50*2/0.339,breaks=20)
abline(v=4.47, lty=2)

Have a quick check:

magimage(psf_model)

Let's try fitting all the sources as Sersic profiles!

group1_modellist=list(
  pointsource=list(
    xcen= group1segstats$xmax[1:3],
    ycen= group1segstats$ymax[1:3],
    mag= group1segstats$mag[1:3]
  ),
  sersic=list(
    xcen= group1segstats$xmax[4:5],
    ycen= group1segstats$ymax[4:5],
    mag= group1segstats$mag[4:5],
    re= group1segstats$R50[4:5]/0.339,
    nser= rep(1, 2),
    ang= group1segstats$ang[4:5],
    axrat= group1segstats$axrat[4:5],
    box= rep(0, 2)
  )
)

group1_interval=list(
  pointsource=list(
    xcen= rep(list(c(-1,1)),3),
    ycen= rep(list(c(-1,1)),3),
    mag= rep(list(c(-2,2)),3)
  ),
  sersic=list(
    xcen= rep(list(c(-1,1)),2),
    ycen= rep(list(c(-1,1)),2),
    mag= rep(list(c(-2,2)),2),
    re= rep(list(c(0.5,2)),2),
    nser= rep(list(c(0.5,10)),2),
    ang= rep(list(c(-180,360)),2),
    axrat= rep(list(c(0.1,1)),2),
    box= rep(list(c(-1,1)),2)
  )
)
for(i in 1:3){
  group1_interval$pointsource$xcen[[i]]=group1_modellist$pointsource$xcen[i]+group1_interval$pointsource$xcen[[i]]
  group1_interval$pointsource$ycen[[i]]=group1_modellist$pointsource$ycen[i]+group1_interval$pointsource$ycen[[i]]
  group1_interval$pointsource$mag[[i]]=group1_modellist$pointsource$mag[i]+group1_interval$pointsource$mag[[i]]
}
for(i in 1:2){
  group1_interval$sersic$xcen[[i]]=group1_modellist$sersic$xcen[i]+group1_interval$sersic$xcen[[i]]
  group1_interval$sersic$ycen[[i]]=group1_modellist$sersic$ycen[i]+group1_interval$sersic$ycen[[i]]
  group1_interval$sersic$mag[[i]]=group1_modellist$sersic$mag[i]+group1_interval$sersic$mag[[i]]
  group1_interval$sersic$re[[i]]=group1_modellist$sersic$re[i]*group1_interval$sersic$re[[i]]
}

group1_tofit=list(
  pointsource=list(
    xcen= rep(TRUE, 3),
    ycen= rep(TRUE, 3),
    mag= rep(TRUE, 3)
  ),
  sersic=list(
    xcen= rep(TRUE, 2),
    ycen= rep(TRUE, 2),
    mag= rep(TRUE, 2),
    re= rep(TRUE, 2),
    nser= rep(TRUE, 2),
    ang= rep(TRUE, 2),
    axrat= rep(TRUE, 2),
    box= rep(FALSE, 2)
  )
)

group1_tolog=list(
  pointsource=list(
    xcen= rep(FALSE, 3),
    ycen= rep(FALSE, 3),
    mag= rep(FALSE, 3)
  ),
  sersic=list(
    xcen= rep(FALSE, 2),
    ycen= rep(FALSE, 2),
    mag= rep(FALSE, 2),
    re= rep(TRUE, 2),
    nser= rep(TRUE, 2),
    ang= rep(FALSE, 2),
    axrat= rep(TRUE, 2),
    box= rep(FALSE, 2)
  )
)

sigma=profoundMakeSigma(image$imDat, objects=profound$objects, gain=0.5, sky=profound$sky, skyRMS =profound$skyRMS, plot=TRUE)

group1_Data=profitSetupData(image$imDat-profound$sky, sigma=sigma, modellist=group1_modellist, tofit=group1_tofit, tolog=group1_tolog, intervals=group1_interval, magzero=30, algo.func='optim', psf=psf_model, region=matrix(profound$segim %in% unlist(profound$group$groupsegID[1,'segID'])[1:5], 356), verbose=FALSE)

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

Next we will run a simple optimisation scheme (because we do not have all day...)

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

And check the results!

profitLikeModel(parm=group1_fit$par, Data=group1_Data, makeplots=TRUE, plotchisq=TRUE)

That is looking pretty pretty pretty good. A bit of residual structure in the central part, but that is as we might expect seeing as we did not attempt to create a more complex (e.g. elliptical and boxy) PSF.

We can compare the inputs and outputs easily too:

group1_modellist_fit=profitRemakeModellist(parm=group1_fit$par, Data=group1_Data)$modellist
group1_modellist_fit

We can also see the change in a nicely formatted way if we use the relist function in R:

relist(unlist(group1_modellist)-unlist(group1_modellist_fit), skeleton = group1_modellist)

The stars have migrated position a small amount, as have the extended sources. The extended sources have both hit a position limit, so it might be sensible to expand the allowed x and y centre intervals. Other parameters are barely changed, e.g. the initial ProFound guesses for the Re (re) and angle (ang) are actually very good. The axial ratios for the extended sources has changed a bit, but this should be expected since these are highly elliptical sources and the convolution with the PSF will be puffing them out when measured without accounting for the PSF in ProFound.

The two bright PSF magnitudes have changed the most, which is not surprising since they were blended together. It is notable that the raw magnitudes calculated by ProFound are actually nearer to the fit values than the rotstats variants. In practice, there is rarely an occasion where the raw segmented magnitudes are not a good enough a starting point for optimisation in ProFit, since it is in practice hard for them to be more than a factor two wrong in flux, which is easily accurate enough for a sensible starting point with ProFit.