ProFit: Found 2 Fit

Aaron Robotham

2021-02-23

In this vignette we will demonstrate fitting in the most automated mode available in ProFit. This is very high level, but hopefully usefully so in many cases. Users who need finer level control should look at some of the other vignettes that break various phases we do automatically here into more manual and explicit steps.

Get the latest version of ProFound and ProFit:

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

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

evalglobal = TRUE

First load the libraries we need:

library(ProFit)
#> Loading required package: FITSio
library(ProFound)
#> Loading required package: Rcpp
library(magicaxis)
library(Highlander)

Load the data:

image = readFITS(system.file("extdata", 'VIKING/mystery_VIKING_Z.fits', package="ProFound"))$imDat

And take a look at what we have got:

magimage(image, hicut=1)

It is Z-band data from the VIKING survey. You can see that there are two very bright stars next to our galaxy of interest. At least the nearer of the two is producing enough extended flux that you profiling might be compromised unless we account for it.

Making an Automatic PSF

If we do not have a good model for the PSF, we can use the high level fairly automated function profitAllStarDoFit to find multiple reasonable stars within the frame (6 by default) and fit them all with a common Moffat model for our PSF. This takes about 7 minutes to run on a modern single CPU.

stars_auto = profitAllStarDoFit(image, magzero=30, SBdilate=2, rough=TRUE)
#> Running Found2Fit
#>     Running ProFound
#> Running Highander
#> Iteration 1
#> CMA 1: -1390.597 25.472 77.021 26.352 76.413 25.722 25.327 76.394 77.211 20.875 21.897 22.154 22.366 0.619 0.391
#> 
#> Laplace's Demon was called on Tue Feb 23 17:52:08 2021
#> 
#> Performing initial checks...
#> WARNING: Matrix mask may be rank-deficient.
#> WARNING: Matrix segim may be rank-deficient.
#> WARNING: Matrix region may be rank-deficient.
#> WARNING: Matrix calcregion may be rank-deficient.
#> Algorithm: Componentwise Hit-And-Run Metropolis 
#> 
#> Laplace's Demon is beginning to update...
#> Iteration: 100,   Proposal: Componentwise,   LP: -1392.7
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.
#> Iteration 2
#> 
#> Laplace's Demon was called on Tue Feb 23 17:52:13 2021
#> 
#> Performing initial checks...
#> WARNING: Matrix mask may be rank-deficient.
#> WARNING: Matrix segim may be rank-deficient.
#> WARNING: Matrix region may be rank-deficient.
#> WARNING: Matrix calcregion may be rank-deficient.
#> Algorithm: Componentwise Hit-And-Run Metropolis 
#> 
#> Laplace's Demon is beginning to update...
#> Iteration: 100,   Proposal: Componentwise,   LP: -1394.3
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.

And check the resulting fit:

profitLikeModel(stars_auto$parm, stars_auto$Data, makeplots = TRUE, plotchisq = TRUE)
#> Warning in dt(x, tdof): NaNs produced

#> $LP
#> [1] -1390.597
#> 
#> $Dev
#> [1] 2781.194
#> 
#> $Monitor
#>           LL           LP tend.elapsed 
#>    -1390.597    -1390.597       12.577 
#> 
#> $yhat
#> [1] 1
#> 
#> $parm
#> moffat.xcen1 moffat.xcen2 moffat.xcen3 moffat.xcen4 moffat.ycen1 moffat.ycen2 
#>   25.4717075   77.0208849   26.3522732   76.4128888   25.7223373   25.3272761 
#> moffat.ycen3 moffat.ycen4  moffat.mag1  moffat.mag2  moffat.mag3  moffat.mag4 
#>   76.3944542   77.2109552   20.8751854   21.8968733   22.1541462   22.3662960 
#> moffat.fwhm1  moffat.con1 
#>    0.6191344    0.3910980

And we can look at the model PSF too:

magimage(stars_auto$psf)

And look at the model parameters too (if you are interested):

stars_auto$psf_modellist
#> $moffat
#> $moffat$xcen
#> [1] 25.5
#> 
#> $moffat$ycen
#> [1] 25.5
#> 
#> $moffat$mag
#> [1] 0
#> 
#> $moffat$fwhm
#> [1] 4.160394
#> 
#> $moffat$con
#> [1] 2.460923
#> 
#> $moffat$ang
#> [1] 0
#> 
#> $moffat$axrat
#> [1] 1

Fit a Galaxy

We can now attempt a fit of our target galaxy near the middle of the image. By default this is what is targeted by profitDoFit, but we will specify the loc explicitly to make the process clear. For the first example we will fit just a single Sersic component (the default). The function automatically detects nearby sources and models them a single Sersic profiles to reduce the contamination in our target source.

example_SS = profitDoFit(image=image, loc=c(178,178), psf=stars_auto$psf, magzero=30,
                         SBdilate=2, Ncomp=1)
#> Running Found2Fit
#>     Running ProFound
#> Running Highander
#> Iteration 1
#> CMA 1: -1714.392 29.183 21.93 18.435 1.119 -0.07 45.25 -0.844
#> 
#> Laplace's Demon was called on Tue Feb 23 17:52:25 2021
#> 
#> Performing initial checks...
#> WARNING: Matrix mask may be rank-deficient.
#> WARNING: Matrix segim may be rank-deficient.
#> WARNING: Matrix psf may be rank-deficient.
#> WARNING: Matrix region may be rank-deficient.
#> WARNING: Matrix calcregion may be rank-deficient.
#> Algorithm: Componentwise Hit-And-Run Metropolis 
#> 
#> Laplace's Demon is beginning to update...
#> Iteration: 100,   Proposal: Componentwise,   LP: -1710.3
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.
#> LD Mode 1: -1706.872 29.076 22.053 18.44 1.124 -0.075 44.396 -0.843
#> Iteration 2
#> CMA 2: -1706.721 29.01 22.125 18.441 1.119 -0.067 44.214 -0.841
#> 
#> Laplace's Demon was called on Tue Feb 23 17:52:37 2021
#> 
#> Performing initial checks...
#> WARNING: Matrix mask may be rank-deficient.
#> WARNING: Matrix segim may be rank-deficient.
#> WARNING: Matrix psf may be rank-deficient.
#> WARNING: Matrix region may be rank-deficient.
#> WARNING: Matrix calcregion may be rank-deficient.
#> Algorithm: Componentwise Hit-And-Run Metropolis 
#> 
#> Laplace's Demon is beginning to update...
#> Iteration: 100,   Proposal: Componentwise,   LP: -1709.2
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.
#> LD Mode 2: -1706.496 29.01 22.125 18.441 1.119 -0.067 44.398 -0.841

And check the resulting fit:

profitLikeModel(example_SS$parm, example_SS$Data, makeplots = TRUE, plotchisq = TRUE)
#> Warning in dt(x, tdof): NaNs produced

#> $LP
#> [1] -1706.496
#> 
#> $Dev
#> [1] 3412.993
#> 
#> $Monitor
#>           LL           LP tend.elapsed 
#>    -1706.496    -1706.496       39.296 
#> 
#> $yhat
#> [1] 1
#> 
#> $parm
#>  sersic.xcen  sersic.ycen   sersic.mag    sersic.re  sersic.nser   sersic.ang 
#>   29.0100415   22.1249149   18.4406096    1.1188187   -0.0665719   44.3980868 
#> sersic.axrat 
#>   -0.8411036

We can try to improve the fit by doing a 2 component fit. By setting Ncomp=1.5 we will fit the bulge as an unresolved PSF:

example_Bpsf = profitDoFit(image=image, loc=c(178,178), psf=stars_auto$psf, magzero=30,
                           SBdilate=2, Ncomp=1.5)
#> Running Found2Fit
#>     Running ProFound
#> Running Highander
#> Iteration 1
#> CMA 1: -1802.344 28.84 22.351 18.296 1.226 44.447 -0.949 22.922
#> 
#> Laplace's Demon was called on Tue Feb 23 17:52:51 2021
#> 
#> Performing initial checks...
#> WARNING: Matrix mask may be rank-deficient.
#> WARNING: Matrix segim may be rank-deficient.
#> WARNING: Matrix psf may be rank-deficient.
#> WARNING: Matrix region may be rank-deficient.
#> WARNING: Matrix calcregion may be rank-deficient.
#> Algorithm: Componentwise Hit-And-Run Metropolis 
#> 
#> Laplace's Demon is beginning to update...
#> Iteration: 100,   Proposal: Componentwise,   LP: -1793.6
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.
#> LD Mode 1: -1788.907 28.885 22.26 18.283 1.221 44.388 -0.96 25.389
#> Iteration 2
#> CMA 2: -1787.691 28.914 22.246 18.283 1.213 44.299 -0.937 25.91
#> 
#> Laplace's Demon was called on Tue Feb 23 17:53:02 2021
#> 
#> Performing initial checks...
#> WARNING: Matrix mask may be rank-deficient.
#> WARNING: Matrix segim may be rank-deficient.
#> WARNING: Matrix psf may be rank-deficient.
#> WARNING: Matrix region may be rank-deficient.
#> WARNING: Matrix calcregion may be rank-deficient.
#> Algorithm: Componentwise Hit-And-Run Metropolis 
#> 
#> Laplace's Demon is beginning to update...
#> Iteration: 100,   Proposal: Componentwise,   LP: -1788.2
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.
#> LD Mode 2: -1787.525 28.967 22.188 18.285 1.212 44.541 -0.926 27.748

And check the resulting fit:

profitLikeModel(example_Bpsf$parm, example_Bpsf$Data, makeplots = TRUE, plotchisq = TRUE)
#> Warning in dt(x, tdof): NaNs produced

#> $LP
#> [1] -1787.525
#> 
#> $Dev
#> [1] 3575.05
#> 
#> $Monitor
#>           LL           LP tend.elapsed 
#>    -1787.525    -1787.525       64.125 
#> 
#> $yhat
#> [1] 1
#> 
#> $parm
#>     sersic.xcen     sersic.ycen      sersic.mag       sersic.re      sersic.ang 
#>      28.9671096      22.1883413      18.2850905       1.2124503      44.5409697 
#>    sersic.axrat pointsource.mag 
#>      -0.9257032      27.7478958

We can try to improve the fit by doing a 2 component fit with a fixed Sersic n=4 bulge.

example_B4 = profitDoFit(image=image, loc=c(178,178), psf=stars_auto$psf, magzero=30, 
                         SBdilate=2, bulge_nser_fit=FALSE, Ncomp=2)
#> Running Found2Fit
#>     Running ProFound
#> Running Highander
#> Iteration 1
#> CMA 1: -1714.906 28.952 22.19 20.179 18.483 1.534 1.134 44.61 -0.914
#> 
#> Laplace's Demon was called on Tue Feb 23 17:53:20 2021
#> 
#> Performing initial checks...
#> WARNING: Matrix mask may be rank-deficient.
#> WARNING: Matrix segim may be rank-deficient.
#> WARNING: Matrix psf may be rank-deficient.
#> WARNING: Matrix region may be rank-deficient.
#> WARNING: Matrix calcregion may be rank-deficient.
#> Algorithm: Componentwise Hit-And-Run Metropolis 
#> 
#> Laplace's Demon is beginning to update...
#> Iteration: 100,   Proposal: Componentwise,   LP: -1710.6
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.
#> LD Mode 1: -1709.231 28.995 22.156 21.665 18.456 0.959 1.134 44.362 -0.879
#> Iteration 2
#> 
#> Laplace's Demon was called on Tue Feb 23 17:53:42 2021
#> 
#> Performing initial checks...
#> WARNING: Matrix mask may be rank-deficient.
#> WARNING: Matrix segim may be rank-deficient.
#> WARNING: Matrix psf may be rank-deficient.
#> WARNING: Matrix region may be rank-deficient.
#> WARNING: Matrix calcregion may be rank-deficient.
#> Algorithm: Componentwise Hit-And-Run Metropolis 
#> 
#> Laplace's Demon is beginning to update...
#> Iteration: 100,   Proposal: Componentwise,   LP: -1713.8
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.

And check the resulting fit:

profitLikeModel(example_B4$parm, example_B4$Data, makeplots = TRUE, plotchisq = TRUE)
#> Warning in dt(x, tdof): NaNs produced

#> $LP
#> [1] -1709.231
#> 
#> $Dev
#> [1] 3418.462
#> 
#> $Monitor
#>           LL           LP tend.elapsed 
#>    -1709.231    -1709.231      107.709 
#> 
#> $yhat
#> [1] 1
#> 
#> $parm
#>  sersic.xcen1  sersic.ycen1   sersic.mag1   sersic.mag2    sersic.re1 
#>    28.9952205    22.1562931    21.6652229    18.4563641     0.9585412 
#>    sersic.re2   sersic.ang2 sersic.axrat2 
#>     1.1338273    44.3618763    -0.8793158