

---
title: "ProFit: PSF Fitting Example"
author: "Aaron Robotham"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{ProFit: PSF Fitting Example}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Prepare the test data

First load the libraries we need:

```{r}
library(knitr)
library(ProFit)
library(FITSio)
```

Next we load a table of data describing GAMA galaxies:

```{r}
data('ExampleInit', package="ProFit")
kable(head(ExampleInit, 10))
```

Here we will use an SDSS example:

```{r}
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:

```{r}
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:

```{r}
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:

```{r, fig.width=5, fig.height=5, eval=F}
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:

```{r, fig.width=5, fig.height=5, eval=F}
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')
```

Visually there is clear evidence of some amount of ellipticity since we can see asymmetry in the vertical cuts made either side of the peak (i.e. the blue lines are to the left of the red lines).

## Setup the fitting data structures

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

```{r}
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:

```{r}
psfmodel=profitMakeModel(modellist,dim=c(25,25))$z
```

```{r, fig.width=5, fig.height=5, eval=F}
magimage(psfmodel)
magimage(abs(psfmodel-psf)/psfmodel,zlim=c(0,1))
```

We will fit everything:

```{r}
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:

```{r}
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:

```{r}
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:

```{r}
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:

```{r,  fig.width=7, fig.height=5, eval=F}
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:

```{r, eval=FALSE, fig.width=5, fig.height=5}
optimfit=optim(Data$init, profitLikeModel, method='BFGS', Data=Data, control=list(fnscale=-1))
```

Check the final result:

```{r, eval=FALSE, fig.width=5, fig.height=5}
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!):

```{r, eval=FALSE}
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:

```{r, eval=FALSE}
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:

```{r, eval=FALSE, fig.width=5, fig.height=5}
PSFmodellist=profitRemakeModelList(optimfit$par, modellist, tofit=tofit, tolog=tolog)
PSFmodellist$moffat$xcen=12.5
PSFmodellist$moffat$xcen=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

```{r}
useID=ExampleIDs[1]
image = readFITS(system.file("extdata", paste(datasource,'/',useID,'fitim.fits',sep=''), package="ProFit"))$imDat
mask = readFITS(system.file("extdata", paste(datasource,'/',useID,'mskim.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):

```{r}
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.

```{r}
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):

```{r, fig.width=5, fig.height=5}
magimage(profitMakeModel(modellist,dim=dim(image)))
```

The original image:

```{r, fig.width=5, fig.height=5}
magimage(image)
```

The convolved model (with PSF):

```{r, fig.width=5, fig.height=5}
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):

```{r}
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:

```{r}
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:

```{r}
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:

```{r}
Data=profitSetupData(image=image, sigma=sigma, segim=segim, mask=mask, 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:

```{r, eval=FALSE}
optimfitMod=optim(Data$init, profitLikeModel, method='BFGS', Data=Data,
                  control=list(fnscale=-1))
```

The best `optim` BFGS fit is given by:

```{r, eval=FALSE}
optimfitMod$par
```

Check it out:

```{r, eval=FALSE}
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'))
```

```{r, eval=FALSE}
modeloptim=profitRemakeModelList(optimfitMod$par,Data$modellist,Data$tofit,Data$tolog)
profitEllipsePlot(Data,modeloptim,pixscale=0.339,FWHM=1,SBlim=26)
```

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

```{r}
Data=profitSetupData(image=image, sigma=sigma, segim=segim, mask=mask, 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:

```{r, eval=FALSE}
optimfitEmp=optim(Data$init, profitLikeModel, method='BFGS', Data=Data,
                  control=list(fnscale=-1))
```

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

```{r, eval=FALSE}
optimfitEmp$par
```

Check it out:

```{r, eval=FALSE, fig.width=5, fig.height=5}
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'))
```

```{r, eval=FALSE, fig.width=5, fig.height=5}
modeloptim=profitRemakeModelList(optimfitEmp$par,Data$modellist,Data$tofit,Data$tolog)
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.