profitSetupData: Setup ProFit Data
In ProFit: Fit Projected 2D Profiles to Galaxy Images

Description Usage Arguments Details Value Notes Author(s) See Also Examples

This is a utility function to get the user inputs in the format required for model optimisation / fitting. It will format the PSF (if supplied) and benchmark the available convolution methods, caching any data required for efficient convolution (such as the PSF FFT). This function does all of the book-keeping required to convert the user data into the format required by ProFit.

profitSetupData(image, region, sigma, segim, mask, modellist, tofit, tolog, priors,
intervals, constraints, psf=NULL, psfdim=dim(psf), finesample=1L, psffinesampled=FALSE,
magzero=0, algo.func='LA', like.func="norm", magmu=FALSE, verbose=FALSE,
omp_threads = NULL, openclenv=NULL, openclenv_int=openclenv, openclenv_conv=openclenv,
nbenchmark=0L, nbenchint=nbenchmark, nbenchconv=nbenchmark,
benchintmethods=c("brute"), benchconvmethods = c("brute","fftw"),
benchprecisions="double", benchconvprecisions=benchprecisions,
benchintprecisions=benchprecisions,
benchopenclenvs = profitGetOpenCLEnvs(make.envs = TRUE),
printbenchmark=FALSE, printbenchint=printbenchmark, printbenchconv=printbenchmark)

`image`	Image matrix; required, the galaxy image we want to fit a model to. The galaxy should be approximately central within this image.
`region`	Logical matrix; optional, specifying the parts of the image to be used for fitting this will override the combination of segim and mask that is used otherwise. If provided this matrix must be the same dimensions as image. Can be integer 1/0 or boolean TRUE/FALSE type logic.
`sigma`	Sigma matrix; optional, the measurement errors per pixel (expressed in terms of sigma). This matrix must be the same dimensions as image.
`segim`	Segmentation matrix; optional, the full segmentation map of the image. If region is not provided then value of the central pixel is used to select the segmented pixels of the galaxy we want to fit. The log-likelihood is then computed using only these pixels. This matrix must be the same dimensions as image.
`mask`	Logical matrix; optional, non galaxy parts of the image to mask out, where 1 means mask out and 0 means use for analysis. If region is not provided then 0 values are used to define the common area to use for fitting. This matrix must be the same dimensions as image.
`modellist`	List; required, the initial model list that describes the analytic model to be created. Can contain an analytical PSF model as well. See Details.
`tofit`	Logical list, optional, using exactly the same list structure as modellist. This flags which parameters of the model list should be fitted. Parameters which are not fitted will inherit their values from modellist. NA values mean the parameter will inherit the value of the previous parameter. In practice this is used to force one parameter (like xcen) to be inherited by multiple Sersic profiles, i.e. we want them to share the same centre, which we fit. The first element of the vector should be TRUE in this case, with the joined profiles set to NA. See Details.
`tolog`	Logical list; optional, using exactly the same list structure as modellist. This flags which parameters of the model list should be fitted in log space (i.e. only relevant to set this if the parameter is being fitted). Parameters like size (re) and axial ratio (axrat) are more naturally explored in log-space, so these should typically be set to true. See Details.
`priors`	Function; optional, that takes the new trial modellist (strictly the first argument) and then the initial modellist (strictly second argument) and returns the log-likelihood of the priors. As long as the returned output if a single summed log-likelihood, there is no restriction on what happens internally to the function. You can also parse additional values to be used internally (say Normal sd, as in the galaxy fitting vignette). This very simple or very complex conditional priors can be specified using R functions. See vignettes for an example. If left empty priors will not be used when computing likelihoods.
`intervals`	List; optional, interval limits for each parameter, using a similar list structure to modellist. The limits should be specified as length 2 vectors: c(low, high) in linear parameter space (no matter if tolog is TRUE for this parameter). See Vignettes and Details.
`constraints`	Function; optional, takes the new trial modellist and returns a list with exactly the same structure. This exists for the purpose of allowing complex relationships between parameters. A simple example is given in the Vignettes of not allowing the bulge Re to become larger than the disk Re. You could also use it to specify the offset of, e.g., a Ferrer profile to be linked to that of the Sersic bulge. Your imagination is the limit, as long as the basic structure returns has the same skeleton as modellist.
`psf`	Matrix; optional. An empirical point spread function (PSF) image matrix that ProFit will use to convolve the image, as an alternative to defining an analytical PSF in modellist. This should have odd sizes in each dimension. If the dimension has an even size then the function will internally interpolate it onto an odd sized grid 1 element larger. `profitSetupData` forces negative values to equal 0. During any convolution `profitConvolvePSF` will force the sum of the pixels to equal 1 to ensure flux conservation during convolution of the model image.
`psfdim`	Numeric; optional. Dimensions of the PSF image to generate when fitting an analytic PSF and convolving extended sources. Defaults to the dimensions of psf. Ignored if there are no extended sources or analytic PSF.
`finesample`	An integer factor to determine how much finer of a grid the model image and PSF should be evaluated on. Because the PSF is discretized, convolution introduces additional discretization of the model, diminishing the accuracy of the convolved model. If this parameter is set to an integer greater than one, the model and PSF (but see psffinesampled) will be upsampled prior to convolution, and then downsampled after convolution. The fine sampling factor must be an integer to avoid non-integral re-binning artefacts when downsampling. Large finesample factors will significantly increase convolution time and accuracy, while moderately increasing model generation time and accuracy, so it is recommended to set nbenchmark to at least a few when using this option.
`psffinesampled`	Logical, is the provided PSF already fine-sampled? If this flag is set and an empirical PSF is provided, it will not be interpolated even if finesample is greater than unity.
`magzero`	The magnitude zero point, where values become scaled by the standard scale=10^(-0.4*(mag-magzero)).
`algo.func`	Character string; the fitting functions being used. Allowed options are "optim", "CMA", "LA" and "LD". `profitLikeModel` uses the value of algo.func in the profit.data object to determine the type of output generated for fitting purposes (see `profitSetupData` for details). If this flag is set to either "optim" or "CMA" then it will output the log-likelihood as a single value. If set to "LA" or "LD" then a more complex structure as expected by `LaplaceApproximation` and `LaplacesDemon` (see details for these functions). In practice the simple log-likelihood scalar output as given by setting to "optim" or "CMA" is useful for a large number of maximisation algorithms available within R. In practice the user must ensure that this option is set correctly for the higher level function used to fit the image data.
`like.func`	Character string specifying the likelihood distribution function to use. Chi-Squared "chisq", Normal "norm" (default), Poisson "pois" and Student-T "t" are the currently supported options. Poisson uses the Cash (or C) statistic, and can be accessed identically using "cash" (or "c"). The choice of the Student-T is probably sensible in the regime where the model is not a perfect reflection of the data- i.e. there are asymmetric or spiral features that the models in ProFit will never be able to reproduce. These can cause high tension when using Normal statistics, but the use of the Student-T (with more mass in the distant wings) reduces the dominance of poorly fitting and un-fitable regions. The degrees of freedom (DoF) for the Student-T are evaluated from the data and model directly so as to maximise the likelihood. If the model is an excellent fit than Normal likelihoods are preferred, and this is the default.
`magmu`	Logical vector. If TRUE then the mag parameter in the input modellist list is interpreted as the mean surface brightness within Re in units of mag/pix^2. If this is of length 1 then all mag values will be interpreted in the same sense, otherwise it should be the same length as the number of components being generated. If FALSE mag is taken to mean total magnitude of the integrated profile. Using this flag might be useful for disk components since they occupy and relatively narrow range in surface brightness, but can have essentially any total magnitude.
`verbose`	Logical; if TRUE then the value of parameters currently being assessed will be printed to screen. Useful for prototyping, but typically this produces a lot of screen output and can slow down the fitting process.
`omp_threads`	An integer indicating the number of threads to use to evaluate radial profiles. If not given only one thread is used. openclenv has precedence over this option, so if both are given then OpenCL evaluation takes place.
`openclenv`	If NULL (default) then the CPU is used to compute the profile. If openclenv is a legal pointer to a graphics card of class externalptr then that card will be used to make a GPU based model. This object can be obtained from the `profitOpenCLEnv` function directly. If openclenv='get' then the OpenCL environment is obtained from running `profitOpenCLEnv` with default values (which are usually reasonable).
`openclenv_int`	The OpenCL environment to use for integrating profiles. Defaults to the value specified in openclenv.
`openclenv_conv`	The OpenCL environment to use for PSF convolution. Defaults to the value specified in openclenv.
`nbenchmark`	Integer; the number of times to benchmark the speed of the available convolution and integration methods. The results of this benchmarking are saved, along with the optimal method.
`nbenchint`	Integer; the number of times to benchmark the speed of the available profile integration methods. The results of this benchmarking are saved, along with the optimal benchmarking method. Defaults to the value specified in nbenchmark.
`nbenchconv`	Integer; the number of times to benchmark the speed of the available convolution methods. The results of this benchmarking are saved, along with the optimal method and any additional data required for efficient convolution (such as the FFT of the PSF, if it is not variable). Defaults to the value specified in nbenchmark.
`benchintmethods`	List of strings specifying which profile integration methods to benchmark. See `profitBenchmark` for details.
`benchconvmethods`	List of strings specifying which convolution methods to benchmark. See `profitBenchmark` for details.
`benchprecisions`	List of floating point precisions to benchmark. Available options are "single" and "double". Defaults to "double", which should be used unless you are certain that single-precision roundoff errors are not important.
`benchintprecisions`	List of floating point precisions to benchmark profile integration with. Available options are "single" and "double". Defaults to benchprecisions.
`benchconvprecisions`	List of floating point precisions to benchmark convolution with. Available options are "single" and "double". Defaults to benchprecisions.
`benchopenclenvs`	List of OpenCL environments to benchmark. Defaults to all available environments. The optimal environment will then be used for openclenvint and openclenvconv, overriding any values set there.
`printbenchmark`	Logical; flag to output a summary of benchmarking results. Default false.
`printbenchint`	Logical; flag to output a summary of profile integration benchmarking results. Defaults to printbenchmark.
`printbenchconv`	Logical; flag to output a summary of convolution benchmarking results. Defaults to printbenchmark.

A legal model list (modellist) has the structure of, e.g., list(sersic, ferrer, psf, sky). At least one of sersic, coresersic, moffat, ferrer, king, pointsource, psf or sky should be present. Each of these is itself a list which contain vectors for each relevant parameter. All these vectors should be the same length for each type of model structure.

The parameters that must be specified for sersic are:

xcen: Vector; x centres of the 2D Sersic profiles (can be fractional pixel positions).
ycen: Vector; y centres of the 2D Sersic profiles (can be fractional pixel positions).
mag: Vector; total magnitudes of the 2D Sersic profiles. Converted to flux using 10^(-0.4*(mag-magzero)).
re: Vector; effective radii of the 2D Sersic profiles
nser: Vector; the Sersic indices of the 2D Sersic profiles
ang: Vector; the orientation of the major axis of the profile in degrees. When plotted as an R image the angle (theta) has the convention that 0= | (vertical), 45= \, 90= - (horizontal), 135= /, 180= | (vertical). Values outside the range 0 <= ang <= 180 are allowed, but these get recomputed as ang = ang %% 180.
axrat: Vector; axial ratios of Sersic profiles defined as minor-axis/major-axis, i.e. 1 is a circle and 0 is a line.
box: Vector; the boxiness of the Sersic profiles that trace contours of iso-flux, defined such that r[mod]=(x^(2+box)+y^(2+box))^(1/(2+box)). When box=0 the iso-flux contours will be normal ellipses, but modifications between -1<box<1 will produce visually boxy distortions. Negative values have a pin-cushion effect, whereas positive values have a barrel effect (the major and minor axes staying fixed in all cases).

The parameters that must be specified for coresersic are:

xcen: Vector; x centres of the 2D Sersic profiles (can be fractional pixel positions).
ycen: Vector; y centres of the 2D Sersic profiles (can be fractional pixel positions).
mag: Vector; total magnitudes of the 2D Sersic profiles. Converted to flux using 10^(-0.4*(mag-magzero)).
re: Vector; effective radii of the 2D Sersic profiles
rb: Vector; transition radius of the Sersic profile (from inner power-law to outer Sersic).
nser: Vector; the Sersic indices of the 2D Sersic profiles
a: Vector; strength of transition from inner core to outer Sersic. Larger +ve means sharper.
b: Vector; the inner power-law of the Core-Sersic. Less than 1 is an increasingly flat core.
ang: Vector; the orientation of the major axis of the profile in degrees. When plotted as an R image the angle (theta) has the convention that 0= | (vertical), 45= \, 90= - (horizontal), 135= /, 180= | (vertical). Values outside the range 0 <= ang <= 180 are allowed, but these get recomputed as ang = ang %% 180.
axrat: Vector; axial ratios of Sersic profiles defined as minor-axis/major-axis, i.e. 1 is a circle and 0 is a line.
box: Vector; the boxiness of the Sersic profiles that trace contours of iso-flux, defined such that r[mod]=(x^(2+box)+y^(2+box))^(1/(2+box)). When box=0 the iso-flux contours will be normal ellipses, but modifications between -1<box<1 will produce visually boxy distortions. Negative values have a pin-cushion effect, whereas positive values have a barrel effect (the major and minor axes staying fixed in all cases).

The parameters that must be specified for moffat are:

xcen: Vector; x centres of the 2D Moffat profiles (can be fractional pixel positions).
ycen: Vector; y centres of the 2D Moffat profiles (can be fractional pixel positions).
mag: Vector; total magnitudes of the 2D Moffat profiles. Converted to flux using 10^(-0.4*(mag-magzero)).
fwhm: Vector; full width half max of the Moffat function.
con: Vector; concentration parameter for Moffat functions. Must be larger than 1.
ang: Vector; the orientation of the major axis of the profile in degrees. When plotted as an R image the angle (theta) has the convention that 0= | (vertical), 45= \, 90= - (horizontal), 135= /, 180= | (vertical). Values outside the range 0 <= ang <= 180 are allowed, but these get recomputed as ang = ang %% 180.
axrat: Vector; axial ratios of Moffat profiles defined as minor-axis/major-axis, i.e. 1 is a circle and 0 is a line.
box: Vector; the boxiness of the Moffat profiles that trace contours of iso-flux, defined such that r[mod]=(x^(2+box)+y^(2+box))^(1/(2+box)). When box=0 the iso-flux contours will be normal ellipses, but modifications between -1<box<1 will produce visually boxy distortions. Negative values have a pin-cushion effect, whereas positive values have a barrel effect (the major and minor axes staying fixed in all cases).

The parameters that must be specified for ferrer or ferrers (either allowed) are:

xcen: Vector; x centres of the 2D Ferrer profiles (can be fractional pixel positions).
ycen: Vector; y centres of the 2D Ferrer profiles (can be fractional pixel positions).
mag: Vector; total magnitudes of the 2D Ferrer profiles. Converted to flux using 10^(-0.4*(mag-magzero)).
rout: Vector; the outer limit of the Ferrer profile. Beyond this radius the profile is evaluated as zero.
a: Vector; the global profile power-law slope. 0 would mean a flat top, and +ve increases in intensity towards the centre.
b: Vector; the strength of the profile truncation as it approaches rout.
ang: Vector; the orientation of the major axis of the profile in degrees. When plotted as an R image the angle (theta) has the convention that 0= | (vertical), 45= \, 90= - (horizontal), 135= /, 180= | (vertical). Values outside the range 0 <= ang <= 180 are allowed, but these get recomputed as ang = ang %% 180.
axrat: Vector; axial ratios of Ferrer profiles defined as minor-axis/major-axis, i.e. 1 is a circle and 0 is a line.
box: Vector; the boxiness of the Ferrer profiles that trace contours of iso-flux, defined such that r[mod]=(x^(2+box)+y^(2+box))^(1/(2+box)). When box=0 the iso-flux contours will be normal ellipses, but modifications between -1<box<1 will produce visually boxy distortions. Negative values have a pin-cushion effect, whereas positive values have a barrel effect (the major and minor axes staying fixed in all cases).

The parameters that must be specified for king are:

xcen: Vector; x centres of the 2D Ferrer profiles (can be fractional pixel positions).
ycen: Vector; y centres of the 2D Ferrer profiles (can be fractional pixel positions).
mag: Vector; total magnitudes of the 2D Ferrer profiles. Converted to flux using 10^(-0.4*(mag-magzero)).
rc: Vector; the core radius of the King profile.
rt: Vector, the truncation radius of the King profile. Beyond this radius the profile is evaluated as zero.
a: Vector; the power-law of the King profile.
ang: Vector; the orientation of the major axis of the profile in degrees. When plotted as an R image the angle (theta) has the convention that 0= | (vertical), 45= \, 90= - (horizontal), 135= /, 180= | (vertical). Values outside the range 0 <= ang <= 180 are allowed, but these get recomputed as ang = ang %% 180.
axrat: Vector; axial ratios of Ferrer profiles defined as minor-axis/major-axis, i.e. 1 is a circle and 0 is a line.
box: Vector; the boxiness of the Ferrer profiles that trace contours of iso-flux, defined such that r[mod]=(x^(2+box)+y^(2+box))^(1/(2+box)). When box=0 the iso-flux contours will be normal ellipses, but modifications between -1<box<1 will produce visually boxy distortions. Negative values have a pin-cushion effect, whereas positive values have a barrel effect (the major and minor axes staying fixed in all cases).

The parameters that must be specified for pointsource (see profitMakePointSource for details) are:

xcen: Vector of x centres of the PSFs (can be fractional pixel positions).
ycen: Vectors of y centres of the PSFs (can be fractional pixel positions).
mag: Vectors of total magnitudes of the PSFs. Converted to flux using 10^(-0.4*(mag-magzero)).

The parameters that may be specified for the psf must be a valid model themselves. Using this option allows users to specify an analytic (e.g. Moffat) PSF.

The parameter that must be specified for sky is:

bg: Value per pixel for the background. This should be the value as measured in the original image, i.e. there is no need to worry about the effect of magzero.

An example of a legal model list structure is:

modellist = list(
sersic = list(
xcen = c(180.5, 50),
ycen = c(90, 50),
mag = c(15, 13),
re = c(140, 50),
nser = c(10, 4),
ang = c(60, 135),
axrat = c(0.5, 0.3),
box = c(2,-2)
),
pointsource = list(
xcen = c(34,10,150),
ycen = c(74,120,130),
mag = c(10,13,16)
),
sky = list(
bg = 3e-12
),
)

The parameters to be fitted are defined in a list with the same format as above:

tofit=list(
sersic=list(
xcen= c(T,NA), #We fit for xcen and tie the two together
ycen= c(T,NA), #We fit for ycen and tie the two together
mag= c(T,T),
#Fit for both re= c(T,T),
#Fit for both nser= c(T,F), #Fit for bulge
ang= c(F,T), #Fit for disk
axrat= c(F,T), #Fit for disk
box= c(F,F)
#Fit for neither ),
pointsource=list(
xcen = c(F,F,F),
ycen = c(F,F,F),
mag = c(F,F,F)
),
sky=list(
bg = F
)
)

Parameters that are better explored in log space are defined in a list with the same format as above:

tolog=list(
sersic=list(
xcen= c(F,F),
ycen= c(F,F),
mag= c(F,F),
re= c(T,T), #re is best fit in log space
nser= c(T,T), #nser is best fit in log space
ang= c(F,F),
axrat= c(T,T), #axrat is best fit in log space
box= c(F,F)
),
psf=list(
xcen = c(F,F,F),
ycen = c(F,F,F),
mag = c(F,F,F)
),
sky=list(
bg = F
)
)

ProFit will only only use the priors function if specified:

priors=function(modellist,modellistinit,sigmas=c(2,2,2,2,5,5,1,1,1,1,30,30,0.3,0.3))
LL=sum(
dnorm(modellist$sersic$xcen,modellist$sersic$xcen,sigmas[1:2],log=TRUE),
dnorm(modellist$sersic$ycen,modellist$sersic$ycen,sigmas[3:4],log=TRUE),
dnorm(modellist$sersic$mag,modellist$sersic$mag,sigmas[5:6],log=TRUE),
dnorm(log10(modellist$sersic$re),log10(modellist$sersic$re),sigmas[7:8],log=TRUE),
dnorm(log10(modellist$sersic$nser),log10(modellist$sersic$nser),sigmas[9:10],log=TRUE),
dnorm(log10(modellist$sersic$axrat),log10(modellist$sersic$axrat),sigmas[13:14],log=TRUE)
)
return=LL

ProFit will only only use the intervals list if specified:

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

ProFit will only only use the constraints function if specified:

constraints=function(modellist)
if(modellist$sersic$re[1]>modellist$sersic$re[2])
modellist$sersic$re[1]=modellist$sersic$re[2]

return=modellist

By ProFit convention the bottom-left part of the bottom-left pixel when plotting the image matrix is c(0,0) and the top-right part of the bottom-left pixel is c(1,1), i.e. the mid-point of pixels are half integer values in x and y.

To confuse things a bit, when R plots an image of a matrix it is transposed and re-ordered vertically to how it appears if you print the matrix directly to screen, i.e. compare print(matrix(1:4,2,2)) and image(matrix(1:4,2,2)). The lowest value (1) is top-left when printed but bottom-left when displayed using image (the red pixel). Both are "correct": the issue is whether you consider the first element of a matrix to be the Cartesian x position (movement in x) or a row element (movement in y). Matrices in maths are always written top-left first where the first argument refers to row number, but images by convention are accessed in a Cartesian sense. Hence [3,4] in a maths matrix means 3 down and 4 right from the top-left, but 3 right and 4 up from the bottom-left in an image.

List; complex structure of class profit.data containing:

`init`	The initial parameters to use for fitting. These are parameters where tofit=TRUE, and are extracted from modellist.
`image`	The specified image matrix.
`mask`	The specified mask matrix.
`sigma`	The specified sigma matrix.
`segim`	The specified segim matrix.
`modellist`	The specified modellist list.
`psf`	The specified psf matrix, if any.
`psftype`	The type of PSF - "analytical" if supplied in modellist, "empirical" if supplied in psf, or "none".
`fitpsf`	Logical flag specifying whether the modellist PSF has any parameters tofit.
`algo.func`	The specified algo.func flag.
`mon.names`	Character vector of parameters to be passed when using the LA/LD algorithms. Defaults to c("LL","LP","dof").
`parm.names`	Character vector of parameter names to be passed when using the LA/LD algorithms.
`N`	The number of pixels that will be used in fitting, i.e. the number of image pixels within the segmentation map, which is the same as sum(region).
`region`	Logical matrix specifying which pixels are inside the fitting region.
`calcregion`	Logical matrix specifying which pixels should have their model values calculated and be convolved by the psf.
`usecalcregion`	Logical specifying whether the calcregion matrix should be used; it may be more efficient not to use it.
`convopt`	List including the optimal convolver object and its OpenCL environment (if any).
`benches`	List containing benchmarking results (if any).
`tofit`	The specified tofit list.
`tolog`	The specified tolog list.
`priors`	The specified priors function.
`intervals`	The specified intervals list.
`constraints`	The specified constraints function.
`like.func`	The specified like.func flag.
`magzero`	The specified magzero scalar.
`finesample`	The specified finesample factor.
`imagedim`	The dimensions of the image matrix.
`verbose`	The specified verbose logical.
`magmu`	The specified magmu logical vector.

One of the list outputs of profitSetupData is the calcregion logical matrix. This tells the model generation and convolution codes whether a particular pixel needs to be considered for fitting purposes. It is computed by convolving the logical region matrix (which itself is the elements of segim containing the galaxy to be fitted) with the psf. Values of the convolved matrix output from profitConvolvePSF above 0 are necessary for accurate likelihood evaluation later, and have their pixel value set to TRUE (or 1). This generally has the visual effect of expanding out the region matrix with a square top-hat kernel the same size as the psf matrix.

Aaron Robotham & Dan Taranu

profitMakeModel, profitSersic, profitCoreSersic, profitMoffat, profitFerrer, profitKing, profitConvolvePSF

# 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 = readFITS(system.file("extdata", paste(datasource,'/',useID,'fitim.fits',sep=''),
package="ProFit"))$imDat
sigma = readFITS(system.file("extdata", paste(datasource,'/',useID,'sigma.fits',sep=''),
package="ProFit"))$imDat
segim = readFITS(system.file("extdata", paste(datasource,'/',useID,'segim.fits',sep=''),
package="ProFit"))$imDat
psf = readFITS(system.file("extdata", paste(datasource,'/',useID,'psfim.fits',sep=''),
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 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 minimal data structure we need for likelihood.

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

# Finally, calcualte the likelihood and make a plot:

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

## Not run: 
# If you're brave and your software/drivers are configured correctly, try benchmarking
# with OpenCL and OpenMP:
openclenvs = profitGetOpenCLEnvs(make.envs = TRUE)

Data=profitSetupData(image=image, sigma=sigma, segim=segim, psf=psf,
modellist=modellist, tofit=tofit, tolog=tolog, magzero=0, algo.func='optim', verbose=TRUE,
nbenchmark = 5L, benchconvmethods = profitAvailableConvolvers(),
benchintmethods = profitAvailableIntegrators(), benchopenclenvs = openclenvs,
printbenchmark = TRUE, omp_threads=4)

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

## End(Not run)