ProFit: Multi-frame and Multi-band Fitting

Aaron Robotham

2021-02-18

Load our useful libraries:

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

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

evalglobal = TRUE

Fitting Simple Offset Images

Here we look at co-fitting images with offset positions.

im1 = Rfits_read_image(system.file("extdata/Offset/Zim1.fits",package="ProFit"), ext=2)$imDat
im2 = Rfits_read_image(system.file("extdata/Offset/Zim2.fits",package="ProFit"), ext=2)$imDat

Here we have two similar Z-band frames taken at slightly offset positions of the same primary source. With ProFit it is possible to fit these simultaneously without resorting to stacking.

magimage(im1)

magimage(im2)

First we should extract the PSFs using the profitAllStarDoFit function (this will take about 5 minutes):

PSF1 = profitAllStarDoFit(im1, magzero=30, rough=TRUE, plot=TRUE)
#> Running Found2Fit
#>     Running ProFound
#> [1] "Summary of used sample:"
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -69.0934   0.4005   1.1222   1.2353   1.9317  24.2812 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 4.102557 1.121245 1.290650 5.289717
#> [1] "Using 2634 out of 2634"
#> Running Highander
#> Iteration 1
#> CMA 1: -6608.818 25.609 76.583 25.351 77.035 25.587 25.94 76.586 76.582 15.194 15.363 17.579 17.525 0.456 0.284
#> 
#> Laplace's Demon was called on Thu Feb 18 16:13:21 2021
#> 
#> Performing initial checks...
#> WARNING: Matrix mask may be rank-deficient.
#> WARNING: Matrix sigma 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: -4424.8
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.
#> LD Mode 1: -4419.846 25.703 76.615 25.97 76.73 25.583 25.945 76.291 76.594 15.177 15.372 17.364 17.404 0.457 0.285
#> Iteration 2
#> 
#> Laplace's Demon was called on Thu Feb 18 16:13:27 2021
#> 
#> Performing initial checks...
#> WARNING: Matrix mask may be rank-deficient.
#> WARNING: Matrix sigma 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: -4421.6
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.
#> LD Mode 2: -4419.756 25.703 76.615 25.97 76.731 25.582 25.946 76.292 76.593 15.178 15.371 17.364 17.405 0.457 0.285

#> [1] "Summary of used sample:"
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> -7.2431 -0.5845  0.1321  0.1095  0.8294  7.6762 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 1.227453 1.049454 1.069986 2.431338
#> [1] "Using 2634 out of 2634"

PSF2 = profitAllStarDoFit(im2, magzero=30, rough=TRUE, plot=TRUE)
#> Running Found2Fit
#>     Running ProFound
#> [1] "Summary of used sample:"
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -59.8426   0.3988   1.1786   1.1924   2.1024  15.2472 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 3.473659 1.244300 1.436738 5.153411
#> [1] "Using 1829 out of 1829"
#> Running Highander
#> Iteration 1
#> CMA 1: -3177.207 25.585 76.279 25.356 76.19 25.463 25.101 76.332 77.048 15.39 17.002 17.158 17.395 0.382 0.273
#> 
#> Laplace's Demon was called on Thu Feb 18 16:13:40 2021
#> 
#> Performing initial checks...
#> WARNING: Matrix mask may be rank-deficient.
#> WARNING: Matrix sigma 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: -3025.5
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.
#> LD Mode 1: -3020.679 25.563 76.289 25.426 76.138 25.443 25.075 76.341 76.977 15.381 16.946 17.185 17.385 0.386 0.272
#> Iteration 2
#> 
#> Laplace's Demon was called on Thu Feb 18 16:13:47 2021
#> 
#> Performing initial checks...
#> WARNING: Matrix mask may be rank-deficient.
#> WARNING: Matrix sigma 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: -3018.5
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.
#> LD Mode 2: -3018.092 25.564 76.299 25.425 76.133 25.453 25.068 76.339 76.984 15.376 16.946 17.183 17.375 0.384 0.27

#> [1] "Summary of used sample:"
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> -7.7484 -0.5688  0.1592  0.1251  0.8445  7.4781 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 1.203135 1.047285 1.064830 2.333613
#> [1] "Using 1829 out of 1829"

We can compare the models:

PSF1$psf_modellist
#> $moffat
#> $moffat$xcen
#> [1] 25.5
#> 
#> $moffat$ycen
#> [1] 25.5
#> 
#> $moffat$mag
#> [1] 0
#> 
#> $moffat$fwhm
#> [1] 2.863622
#> 
#> $moffat$con
#> [1] 1.925433
#> 
#> $moffat$ang
#> [1] 0
#> 
#> $moffat$axrat
#> [1] 1
PSF2$psf_modellist
#> $moffat
#> $moffat$xcen
#> [1] 25.5
#> 
#> $moffat$ycen
#> [1] 25.5
#> 
#> $moffat$mag
#> [1] 0
#> 
#> $moffat$fwhm
#> [1] 2.423221
#> 
#> $moffat$con
#> [1] 1.861939
#> 
#> $moffat$ang
#> [1] 0
#> 
#> $moffat$axrat
#> [1] 1

The seeing in our second image is clearly quite a bit better (FWHM 2.4 versus 2.9). This is also why we should get a better solution by fitting both images simultaneously rather than stacking.

To get most of the way to the correct Datalist structure we need we can use profitFound2Fit:

F2F1 = profitFound2Fit(im1, loc=c(450,500), psf=PSF1$psf, cutbox=c(300,300), Ncomp=2,
                       magzero=30, SBdilate=2, tightcrop=FALSE)
#> No input segim that matches the input image- will create one using ProFound!
#>     Running ProFound
F2F2 = profitFound2Fit(im2, loc=c(450,500), psf=PSF2$psf, cutbox=c(300,300), Ncomp=2,
                       magzero=30, SBdilate=2, tightcrop=FALSE)
#> No input segim that matches the input image- will create one using ProFound!
#>     Running ProFound

The main adjustment we need to make it to tell ProFit what the offset is between the frames and put this into the Data object contained. From the original WCSs the offset from im1 to im2 is know to be [-23.112,-9.697] pixels:

F2F2$Data$offset = c(-23.112,-9.697)

We can now check this using the profitLikeModel function, where the same initial parameters get shifted for the second Data object:

profitLikeModel(F2F1$Data$init, Data=F2F1$Data, makeplots=TRUE)
#> [1] "Summary of used sample:"
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> -6.6857 -0.4531  0.3027  0.4301  1.1260 24.4117 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 1.582811 1.173234 1.237891 2.808601
#> [1] "Using 22881 out of 22881"

#> $LP
#> [1] -51802.61
#> 
#> $Dev
#> [1] 103605.2
#> 
#> $Monitor
#>           LL           LP tend.elapsed 
#>   -51802.607   -51802.607       42.807 
#> 
#> $yhat
#> [1] 1
#> 
#> $parm
#>  sersic.xcen1  sersic.ycen1   sersic.mag1   sersic.mag2    sersic.re1 
#>   161.0384678   153.8342988    14.0231202    14.0231202     1.2186830 
#>    sersic.re2   sersic.ang2 sersic.axrat2 
#>     1.6958043   114.0776665    -0.2732623
profitLikeModel(F2F1$Data$init, Data=F2F2$Data, makeplots=TRUE)
#> [1] "Summary of used sample:"
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> -9.7233 -0.2640  0.4846  0.5947  1.3009 26.2145 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 1.602233 1.154734 1.203201 2.813663
#> [1] "Using 17968 out of 17968"

#> $LP
#> [1] -42751.26
#> 
#> $Dev
#> [1] 85502.52
#> 
#> $Monitor
#>           LL           LP tend.elapsed 
#>   -42751.261   -42751.261       43.033 
#> 
#> $yhat
#> [1] 1
#> 
#> $parm
#>  sersic.xcen1  sersic.ycen1   sersic.mag1   sersic.mag2    sersic.re1 
#>   137.9264678   144.1372988    14.0231202    14.0231202     1.2186830 
#>    sersic.re2   sersic.ang2 sersic.axrat2 
#>     1.6958043   114.0776665    -0.2732623

The solutions look pretty similar now! To simplify things we want to create a Datalist structure that contains the two lists:

Datalist = c(list(F2F1$Data), list(F2F2$Data))

You can now run profitLikeModel more directly:

profitLikeModel(F2F1$Data$init, Data=Datalist, makeplots=TRUE)
#> [1] "Summary of used sample:"
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> -6.6857 -0.4531  0.3027  0.4301  1.1260 24.4117 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 1.582811 1.173234 1.237891 2.808601
#> [1] "Using 22881 out of 22881"