# 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]
box = c(160,160)
image = magcutout(readFITS(system.file("extdata", paste(datasource,'/',useID,'fitim.fits',sep=''),
package="ProFit"))$imDat, box = box)$image
sigma = magcutout(readFITS(system.file("extdata", paste(datasource,'/',useID,'sigma.fits',sep=''),
package="ProFit"))$imDat, box = box)$image
segim = magcutout(readFITS(system.file("extdata", paste(datasource,'/',useID,'segim.fits',sep=''),
package="ProFit"))$imDat, box = box)$image
psf = magcutout(readFITS(system.file("extdata", paste(datasource,'/',useID,'psfim.fits',sep=''),
package="ProFit"))$imDat, box = c(21,21))$image
# 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 together
ycen= c(TRUE,NA), #We fit for ycen and tie the two together
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)
)
)
# Setup the profit.data
openclenvs = data.frame()
Data=profitSetupData(image=image, sigma=sigma, segim=segim, psf=psf,
modellist=modellist, tofit=tofit, tolog=tolog, magzero=0, algo.func='optim', verbose=TRUE,
nbenchmark = 1L, benchconvmethods = "brute",
benchintmethods = "brute", benchopenclenvs = openclenvs,
finesample=4L, printbenchmark = TRUE)
system.time(profitLikeModel(parm=Data$init, Data=Data))
benchmarks = profitDataBenchmark(modellist = Data$modellist, calcregion = Data$calcregion,
imgdim = dim(Data$image), finesample = Data$finesample, psf = Data$psf, fitpsf = Data$fitpsf,
nbenchmark = 1L, benchconvmethods = profitAvailableConvolvers(),
benchintmethods = profitAvailableIntegrators(), benchopenclenvs = openclenvs,
printbenchmark = TRUE)
Data = profitDataSetOptionsFromBenchmarks(Data, benchmarks)
system.time(profitLikeModel(parm=Data$init, Data=Data))
# }
Run the code above in your browser using DataLab