ProFit: Offset Images

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, plot=TRUE)
#> Running Found2Fit
#>     Running ProFound
#> [1] "Summary of used sample:"
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -54.0382   0.3868   1.0965   1.1857   1.9001  22.9720 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 3.801279 1.115354 1.276093 5.403320
#> [1] "Using 2634 out of 2634"
#> Running Highander
#> Iteration 1
#> CMA 1: -6193.77 25.678 76.577 26.219 76.354 25.602 25.943 76.763 76.618 15.185 15.398 17.557 17.601 0.443 0.285
#> 
#> Laplace's Demon was called on Thu Feb  4 15:11:14 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: -4444.5
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.
#> LD Mode 1: -4434.757 25.701 76.614 25.973 76.724 25.581 25.946 76.287 76.604 15.181 15.372 17.372 17.398 0.447 0.282
#> Iteration 2
#> 
#> Laplace's Demon was called on Thu Feb  4 15:12:48 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: -4437.7
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.
#> LD Mode 2: -4431.386 25.696 76.614 25.98 76.728 25.576 25.951 76.287 76.619 15.177 15.368 17.368 17.397 0.446 0.28

#> [1] "Summary of used sample:"
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> -7.6589 -0.5816  0.1373  0.1075  0.8405  6.7160 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 1.231222 1.053084 1.071503 2.452146
#> [1] "Using 2634 out of 2634"

PSF2 = profitAllStarDoFit(im2, magzero=30, plot=TRUE)
#> Running Found2Fit
#>     Running ProFound
#> [1] "Summary of used sample:"
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -39.6172   0.3709   1.1509   1.1267   2.0631  14.6307 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 3.145290 1.230457 1.388358 5.352191
#> [1] "Using 1829 out of 1829"
#> Running Highander
#> Iteration 1
#> CMA 1: -3788.664 25.586 76.287 25.472 76.25 25.397 25.064 76.224 76.739 15.4 16.899 17.193 17.296 0.365 0.266
#> 
#> Laplace's Demon was called on Thu Feb  4 15:14:26 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: -2987.5
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.
#> LD Mode 1: -2987.546 25.564 76.29 25.435 76.146 25.443 25.079 76.358 76.982 15.369 16.934 17.17 17.366 0.359 0.257
#> Iteration 2
#> 
#> Laplace's Demon was called on Thu Feb  4 15:16:18 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: -2984.1
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.
#> LD Mode 2: -2983.492 25.566 76.274 25.418 76.151 25.452 25.079 76.343 76.977 15.369 16.938 17.169 17.374 0.36 0.255

#> [1] "Summary of used sample:"
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -7.29631 -0.62025  0.10776  0.08149  0.80825  6.50246 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 1.191087 1.063010 1.086887 2.372173
#> [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.791671
#> 
#> $moffat$con
#> [1] 1.904349
#> 
#> $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.291604
#> 
#> $moffat$con
#> [1] 1.79823
#> 
#> $moffat$ang
#> [1] 0
#> 
#> $moffat$axrat
#> [1] 1

The seeing in our second image is clearly quite a bit better (FWHM 2.3 versus 2.8). 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)
#>     Running ProFound
F2F2 = profitFound2Fit(im2, loc=c(450,500), psf=PSF2$psf, cutbox=c(300,300), Ncomp=2,
                       magzero=30, SBdilate=2, tightcrop=FALSE)
#>     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. 
#> -24.0255  -0.4638   0.3144   0.4133   1.1256  34.2776 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 1.954969 1.176549 1.233549 3.011684
#> [1] "Using 22881 out of 22881"

profitLikeModel(F2F1$Data$init, Data=F2F2$Data, makeplots=TRUE)
#> [1] "Summary of used sample:"
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -28.9496  -0.2437   0.4894   0.5815   1.3033  36.7655 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 2.028850 1.144899 1.189113 3.014019
#> [1] "Using 17968 out of 17968"

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. 
#> -24.0255  -0.4638   0.3144   0.4133   1.1256  34.2776 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 1.954969 1.176549 1.233549 3.011684
#> [1] "Using 22881 out of 22881"

#> [1] "Summary of used sample:"
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -28.9496  -0.2437   0.4894   0.5815   1.3033  36.7655 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 2.028850 1.144899 1.189113 3.014019
#> [1] "Using 17968 out of 17968"

#> $LP
#> [1] -123231
#> 
#> $Dev
#> [1] 246461.9
#> 
#> $Monitor
#>           LL           LP tend.elapsed 
#>  -123230.953  -123230.953      839.128 
#> 
#> $yhat
#> [1] 1
#> 
#> $parm
#>  sersic.xcen1  sersic.ycen1   sersic.mag1   sersic.mag2    sersic.re1 
#>   162.0384678   154.8342988    14.0231202    14.0231202     1.2186830 
#>    sersic.re2  sersic.nser1  sersic.nser2   sersic.ang2 sersic.axrat2 
#>     1.6958043     0.6020600     0.0000000   114.0776665    -0.2732623

The final thing we need to do is add a couple of elements to the Datalist to make it compatible with LaplacesDemon:

Datalist$mon.names = F2F1$Data$mon.names
Datalist$parm.names = F2F1$Data$parm.names
Datalist$N = F2F1$Data$N

Now we can finally run Highlander to find the best common solution for the two frames (this will take about 5 minutes):

HighFit = Highlander(F2F1$Data$init, Data=Datalist, likefunc=profitLikeModel, ablim=0.5)
#> Iteration 1
#> CMA 1: -69848.951 160.051 153.66 14.018 13.878 1.194 1.666 1.045 -0.221 112.695 -0.334
#> 
#> Laplace's Demon was called on Thu Feb  4 15:18:56 2021
#> 
#> Performing initial checks...
#> Algorithm: Componentwise Hit-And-Run Metropolis 
#> 
#> Laplace's Demon is beginning to update...
#> Iteration: 100,   Proposal: Componentwise,   LP: -68216.2
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.
#> LD Mode 1: -68212.212 159.869 153.476 14.066 13.821 1.197 1.664 1.252 -0.119 113.594 -0.326
#> Iteration 2
#> 
#> Laplace's Demon was called on Thu Feb  4 15:22:19 2021
#> 
#> Performing initial checks...
#> Algorithm: Componentwise Hit-And-Run Metropolis 
#> 
#> Laplace's Demon is beginning to update...
#> Iteration: 100,   Proposal: Componentwise,   LP: -68184.9
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.
#> LD Mode 2: -68184.885 159.862 153.477 14.08 13.819 1.187 1.665 1.704 -0.123 113.592 -0.326

And now we can check the results:

profitLikeModel(HighFit$parm, Data=Datalist, makeplots=TRUE)
#> [1] "Summary of used sample:"
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -5.90366 -0.61882  0.07365  0.05201  0.70091  8.79113 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 1.2060390 0.9737905 0.9882521 2.2523414
#> [1] "Using 22881 out of 22881"

#> [1] "Summary of used sample:"
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -7.38031 -0.55313  0.11253  0.07543  0.71825  8.63349 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 1.2450283 0.9357179 0.9642497 2.3605208
#> [1] "Using 17968 out of 17968"

#> $LP
#> [1] -68184.89
#> 
#> $Dev
#> [1] 136369.8
#> 
#> $Monitor
#>           LL           LP tend.elapsed 
#>   -68184.885   -68184.885     1656.091 
#> 
#> $yhat
#> [1] 1
#> 
#> $parm
#>  sersic.xcen1  sersic.ycen1   sersic.mag1   sersic.mag2    sersic.re1 
#>   159.8616538   153.4772459    14.0802053    13.8186502     1.1871149 
#>    sersic.re2  sersic.nser1  sersic.nser2   sersic.ang2 sersic.axrat2 
#>     1.6652059     0.7242759    -0.1231516   113.5921656    -0.3261178

A similar methodology will work on any number of images, as long as they fit within memory. They should strictly all be the same band and pixel scale however, since currently (as of Jan 2021) ProFit does not allow the main model parameters to vary between images (other than the [X,Y] offset as shown above). Longer term the plan is to allow such functionality.

Free, Fixed, and Follow Fits

It is also possible, but more complicated, to fit some but not all parameters. For this mode the users will need to be keenly aware of which parameters are which within the parm stuctures of Data since they will need to construct their own, and tag which parts should correspond to the relevant model list. If that sounds complicated, then its because it is! With great flexibility comes great complexity…

Imagine we want to do a fit much like the above, but where we want to the two disk Re argument to be free to fit differently in image 2. In this case the image is of the same object in the same band, so there is no real reason to allow the disk Re to fit freely (they should be linked together) but we can treat it as a placeholder for an observation in a different band. In this case we might expect many parameters to be precisely linked ([X,Y] and position angle), but the size (Re) to be different because we expect colour gradients in (at least some) disks. As a simple example we can look at the above Datalist.

Datalist[[1]]$init
#>  sersic.xcen1  sersic.ycen1   sersic.mag1   sersic.mag2    sersic.re1 
#>   162.0384678   154.8342988    14.0231202    14.0231202     1.2186830 
#>    sersic.re2  sersic.nser1  sersic.nser2   sersic.ang2 sersic.axrat2 
#>     1.6958043     0.6020600     0.0000000   114.0776665    -0.2732623
Datalist[[2]]$init
#>  sersic.xcen1  sersic.ycen1   sersic.mag1   sersic.mag2    sersic.re1 
#>   138.9622250   145.2416319    14.0558577    14.0558577     1.2128408 
#>    sersic.re2  sersic.nser1  sersic.nser2   sersic.ang2 sersic.axrat2 
#>     1.6899620     0.6020600     0.0000000   116.3105014    -0.2890468

We now need to make a new init structure that adds on an additional re term for the second image disk (re2b) and a new parmuse structure that tells ProFit where to look for parameters when fittin each image.

Datalist2 = c(list(F2F1$Data), list(F2F2$Data)) #Note we are still using the offset added to F2F2$Data

init_new = F2F1$Data$init
init_new = c(init_new, sersic.re2b = as.numeric(F2F2$Data$init['sersic.re2']))
init_new
#>  sersic.xcen1  sersic.ycen1   sersic.mag1   sersic.mag2    sersic.re1 
#>   162.0384678   154.8342988    14.0231202    14.0231202     1.2186830 
#>    sersic.re2  sersic.nser1  sersic.nser2   sersic.ang2 sersic.axrat2 
#>     1.6958043     0.6020600     0.0000000   114.0776665    -0.2732623 
#>   sersic.re2b 
#>     1.6899620

Datalist2$mon.names = F2F1$Data$mon.names
Datalist2$parm.names = names(init_new) #Careful: This needs to reflect the new longer parm init.
Datalist2$N = F2F1$Data$N

Datalist2[[1]]$init = init_new
Datalist2[[2]]$init = init_new

Datalist2[[1]]$parmuse = 1:10 #We use the first arguments for image 1 fitting
Datalist2[[2]]$parmuse = c(1,2,3,4,5,11,7,8,9,10)

We can sanity check this by passing through profitLikeModel as before:

profitLikeModel(Datalist2[[1]]$init, Data=Datalist2, makeplots=TRUE)
#> [1] "Summary of used sample:"
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -24.0255  -0.4638   0.3144   0.4133   1.1256  34.2776 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 1.954969 1.176549 1.233549 3.011684
#> [1] "Using 22881 out of 22881"

#> [1] "Summary of used sample:"
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -29.0521  -0.2355   0.4975   0.5780   1.2983  36.7003 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 2.023474 1.136721 1.178705 2.994912
#> [1] "Using 17968 out of 17968"

#> $LP
#> [1] -122999.1
#> 
#> $Dev
#> [1] 245998.2
#> 
#> $Monitor
#>           LL           LP tend.elapsed 
#>  -122999.080  -122999.080     1657.423 
#> 
#> $yhat
#> [1] 1
#> 
#> $parm
#>  sersic.xcen1  sersic.ycen1   sersic.mag1   sersic.mag2    sersic.re1 
#>   162.0384678   154.8342988    14.0231202    14.0231202     1.2186830 
#>    sersic.re2  sersic.nser1  sersic.nser2   sersic.ang2 sersic.axrat2 
#>     1.6958043     0.6020600     0.0000000   114.0776665    -0.2732623 
#>   sersic.re2b 
#>     1.6899620

Now we can finally run Highlander to find the best common solution for the two frames (this will take about 5 minutes) where the disk Re are allowed to differ:

HighFit2 = Highlander(Datalist2[[1]]$init, Data=Datalist2, likefunc=profitLikeModel, ablim=0.5)
#> Iteration 1
#> CMA 1: -70114.085 160.053 153.497 13.784 14.038 1.329 1.675 0.935 -0.195 114.138 -0.334 1.665
#> 
#> Laplace's Demon was called on Thu Feb  4 15:25:54 2021
#> 
#> Performing initial checks...
#> Algorithm: Componentwise Hit-And-Run Metropolis 
#> 
#> Laplace's Demon is beginning to update...
#> Iteration: 100,   Proposal: Componentwise,   LP: -69078.2
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.
#> LD Mode 1: -69038.663 159.877 153.486 13.829 13.971 1.311 1.672 1.16 -0.187 113.969 -0.349 1.675
#> Iteration 2
#> 
#> Laplace's Demon was called on Thu Feb  4 15:29:32 2021
#> 
#> Performing initial checks...
#> Algorithm: Componentwise Hit-And-Run Metropolis 
#> 
#> Laplace's Demon is beginning to update...
#> Iteration: 100,   Proposal: Componentwise,   LP: -68895.6
#> 
#> Assessing Stationarity
#> Assessing Thinning and ESS
#> Creating Summaries
#> Creating Output
#> 
#> Laplace's Demon has finished.
#> LD Mode 2: -68895.645 159.861 153.471 13.863 13.956 1.3 1.672 1.387 -0.173 113.917 -0.346 1.675

We still a get a pretty good fit:

profitLikeModel(HighFit2$parm, Data=Datalist2, makeplots=TRUE)
#> [1] "Summary of used sample:"
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -6.36246 -0.61525  0.08388  0.05631  0.71657  8.81574 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 1.2189886 0.9817175 0.9922362 2.2619830
#> [1] "Using 22881 out of 22881"

#> [1] "Summary of used sample:"
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -6.75605 -0.55800  0.12171  0.07642  0.73398  8.69899 
#> [1] "sd / MAD / 1-sig / 2-sig range:"
#> [1] 1.2603495 0.9497893 0.9784550 2.3966272
#> [1] "Using 17968 out of 17968"

#> $LP
#> [1] -68895.64
#> 
#> $Dev
#> [1] 137791.3
#> 
#> $Monitor
#>           LL           LP tend.elapsed 
#>   -68895.645   -68895.645     2540.354 
#> 
#> $yhat
#> [1] 1
#> 
#> $parm
#>  sersic.xcen1  sersic.ycen1   sersic.mag1   sersic.mag2    sersic.re1 
#>   159.8611583   153.4705585    13.8628368    13.9564352     1.2999016 
#>    sersic.re2  sersic.nser1  sersic.nser2   sersic.ang2 sersic.axrat2 
#>     1.6723315     0.7242759    -0.1728549   113.9167555    -0.3464978 
#>   sersic.re2b 
#>     1.6746246

Now we can take a look at the fit parameters. Note that even though they were free to be fitted separately, sersic.re2 and sersic.re2b are very similar!

HighFit2$parm
#>  sersic.xcen1  sersic.ycen1   sersic.mag1   sersic.mag2    sersic.re1 
#>   159.8611583   153.4705585    13.8628368    13.9564352     1.2999016 
#>    sersic.re2  sersic.nser1  sersic.nser2   sersic.ang2 sersic.axrat2 
#>     1.6723315     1.3865879    -0.1728549   113.9167555    -0.3464978 
#>   sersic.re2b 
#>     1.6746246

ProFit: Offset Images

Aaron Robotham

2021-02-04

Fitting Simple Offset Images

Free, Fixed, and Follow Fits