ProFound-Complex-Fit.Rmd

---
title: "ProFound/ProFit: A Complex Fit"
author: "Aaron Robotham"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{ProFound/ProFit: A Complex Fit}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

Get the latest version of **ProFound** and **ProFit**:

```{r, eval=FALSE}
library(devtools)
install_github('asgr/ProFound')
install_github('ICRAR/ProFit')
```

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

```{r}
evalglobal=FALSE
```

First load the libraries we need:

```{r, eval=evalglobal}
library(ProFound)
library(ProFit)
library(Rfits)
library(Rwcs)
```

```{r, eval=evalglobal, fig.width=5, fig.height=5, dpi=40}
image = Rfits_read_image(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**:

```{r, eval=evalglobal}
profound$group$groupsegID[1,]
```

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

```{r, eval=evalglobal, fig.width=5, fig.height=5, dpi=40}
magimage(profound$group$groupim==1)
```

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

```{r, eval=evalglobal}
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.

```{r, eval=evalglobal}
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**.

```{r, eval=evalglobal, fig.width=5, fig.height=5, dpi=40}
maghist(profound$segstats$R50*2/0.339,breaks=20)
abline(v=4.47, lty=2)
```

Have a quick check:

```{r, eval=evalglobal, fig.width=5, fig.height=5, dpi=40}
magimage(psf_model)
```

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

```{r, eval=evalglobal}
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],
    ycen= group1segstats$ymax[4],
    mag= group1segstats$mag[4],
    re= group1segstats$R50[4]/0.339,
    nser= 1,
    ang= group1segstats$ang[4],
    axrat= group1segstats$axrat[4],
    box= 0
  )
)

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= list(c(-1,1)),
    ycen= list(c(-1,1)),
    mag= list(c(-2,2)),
    re= list(c(0.5,2)),
    nser= list(c(0.5,10)),
    ang= list(c(-180,360)),
    axrat= list(c(0.1,1)),
    box= list(c(-1,1))
  )
)
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){
  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= TRUE,
    ycen= TRUE,
    mag= TRUE,
    re= TRUE,
    nser= TRUE,
    ang= TRUE,
    axrat= TRUE,
    box= FALSE
  )
)

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


```{r, eval=evalglobal, fig.width=5, fig.height=5, dpi=40}
sigma=profoundMakeSigma(image$imDat, objects=profound$objects, gain=0.5, sky=profound$sky, skyRMS=profound$skyRMS, plot=TRUE)
```

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

```{r, eval=evalglobal, fig.width=8, fig.height=5, dpi=40}
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...)

```{r, echo=TRUE, message=FALSE, warning=FALSE, eval=evalglobal}
group1_fit=optim(group1_Data$init, profitLikeModel, method='BFGS', Data=group1_Data, control=list(fnscale=-1))
```

And check the results!

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

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

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