ProFit-package: Fit Projected 2D Profiles to Galaxy Images

In ICRAR/ProFit: Fit Projected 2D Profiles to Galaxy Images

Fit Projected 2D Profiles to Galaxy Images

Description

Get data / Define model / ??? / Profit! 'ProFit' is a Bayesian galaxy fitting tool that uses a fast 'C++' image generation library and a flexible interface to a large number of likelihood samplers.

Details

Package:	ProFit
Type:	Package
Version:	2.5.0
Date:	2023-03-16
License:	LGPL-3
Depends:	R (>= 3.0), Rfits (>= 1.8.0), magicaxis (>= 2.0.3)
Imports:	cubature, RColorBrewer, LaplacesDemon, methods, celestial (>= 1.4.1), checkmate
Suggests:	fftw, knitr, rmarkdown, ProFound (>= 1.15.0), sn, Highlander (>= 0.1.7), ProSpect

Author(s)

NA

Maintainer: Aaron Robotham <[e-mail]>

References

Robotham A.S.G., et al., 2017, MNRAS, 466, 1513

Examples

modellist = list(
  sersic = list(
    xcen   = c(180, 60),
    ycen   = c(90, 10),
    mag = c(15, 13),
    re  = c(14, 5),
    nser  = c(3, 10),
    ang  = c(46, 80),
    axrat  = c(0.4, 0.6),
    box = c(0.5,-0.5)
  ),
  pointsource = list(
    xcen   = c(34,10,150),
    ycen   = c(74,120,130),
    mag = c(10,13,16)
  ),
  sky = list(
    bg = 3e-12
  )
)

# Without a PSF provided only the extended sources are shown, with no convolution:

magimage(profitMakeModel(modellist=modellist, dim=c(200,200)))

# With a PSF provided the PSFs are displayed and the extended sources are convolved with
# the PSF:

magimage(profitMakeModel(modellist=modellist, psf=profitMakePointSource(), dim=c(200,200)))

############### Full L-BFGS-B fit example ##############

## Not run: 

# Load ProFit example data

# There are 2 data source options: KiDS or SDSS (the galaxies are the same)

datasource='KiDS' 

# 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:

data('ExampleInit')
ExampleFiles = list.files(system.file("extdata",datasource,package="ProFit"))
ExampleIDs = unlist(strsplit(ExampleFiles[grep('fitim',ExampleFiles)],'fitim.fits'))
print(ExampleIDs)

# There are 10 example galaxies included. Here we run example 1:

useID=ExampleIDs[1]

image = Rfits_read_image(system.file("extdata", paste0(datasource,'/',useID,'fitim.fits'),
package="ProFit"))$imDat
sigma = Rfits_read_image(system.file("extdata", paste0(datasource,'/',useID,'sigma.fits'),
package="ProFit"))$imDat
segim = Rfits_read_image(system.file("extdata", paste0(datasource,'/',useID,'segim.fits'),
package="ProFit"))$imDat
psf = Rfits_read_image(system.file("extdata", paste0(datasource,'/',useID,'psfim.fits'),
package="ProFit"))$imDat

# Very rough model (not meant to look too good yet):

useIDnum=as.integer(strsplit(useID,'G')[[1]][2])
useloc=which(ExampleInit$CATAID==useIDnum)

# For our initial model we treat component 1 as the putitive bulge and componet 2 as
# the putitive 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.

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)
  )
)

# The pure model (no PSF):
magimage(profitMakeModel(modellist,dim=dim(image)))

# The original image:
magimage(image)

# The convolved model (with PSF):
magimage(profitMakeModel(modellist,dim=dim(image),psf=psf))

# What should we be fitting:

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
  )
)

# What parameters should be fitted in log space:

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 interval limits to use when fitting. This is not strictly required, but without
# this we cannot ensure the sampler does not enter unallowed values like negative sizes,
# Sersic indices and axial ratios etc:

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 minimal data structure we need for optimisation. See vignettes for
# more complex examples using priors, and constraints:

Data=profitSetupData(image=image, sigma=sigma, segim=segim,psf=psf,
modellist=modellist, tofit=tofit, tolog=tolog, intervals=intervals, magzero=0,
algo.func='optim', verbose=TRUE)

# This produces a fairly complex R object, but with all the bits we need for fitting,
# e.g. (notice the tolog parameteres are now logged):

Data$init

# These are the parameters we wish to fit for, and we take the initial guesses from the
# model list we provided before.

# We can test how things currently look (we get an output because we set verbose=TRUE
# earlier):

profitLikeModel(parm=Data$init, Data=Data, makeplots=TRUE)

# Let us try optim BFGS:

optimfit=optim(Data$init, profitLikeModel, method='BFGS', Data=Data,
control=list(fnscale=-1))

# The best optim BFGS fit is given by:

optimfit$par

# Check it out:

profitLikeModel(optimfit$par,Data,makeplots=TRUE,whichcomponents=list(sersic=1))
profitLikeModel(optimfit$par,Data,makeplots=TRUE,whichcomponents=list(sersic=2))
profitLikeModel(optimfit$par,Data,makeplots=TRUE,whichcomponents=list(sersic='all'))

modeloptim=profitRemakeModellist(optimfit$par,Data$modellist,Data$tofit,Data$tolog)$modellist
profitEllipsePlot(Data,modeloptim,pixscale=0.2,FWHM=0.5,SBlim=26)

## End(Not run)

ICRAR/ProFit documentation built on Feb. 1, 2024, 9:34 a.m.