GIM2D/GIMFIT2D

Parameters

     in_image = "ngal1"         Name of input image
    out_image = "junk_out"      Name of output image
   (sym_image = "junk_sym")     Name of symmetrized image
   (res_image = "junk_res")     Name of residual image
   (msk_image = "")             Name of pixel mask image
   (sig_image = "none")         Name of sigma image
     (logfile = "junk.log")     Name of log file
   (mdparfile = "gal.mdpar")    Name of parameter value file
   (initparam = yes)            Determine parameters from input image moments (yes/no)?
       (dosym = no)             Symmetrize input image (yes/no)?
       (dobkg = yes)            Recompute background parameters (yes/no)?
     (imscale = 0.1)            Image scale (arcsecond/pixel)
     (psftype = "tinytim")      PSF type (gaussian|tiny_wfpc|tiny_nic|tiny_stis|user)
      (in_psf = "/net/shamash/dsk/extra/simard/bd/test/inimpsf") Name of tinytim/user PSF image
       (c_abs = 0.)             Disk internal absorption coefficient (0-1)
      (seeing = 0.06)           Seeing FWHM (arcseconds)
      (bkgsig = 0.65)           Background sigma (ADU)
       (nsamp = 300)            Number of parameter space samples
     (ccdgain = 7.)             CCD Gain (electrons/DU)
      (nframe = 1)              Number of frames combined
     (metseed = -22029)         Metropolis seed
       (osamp = 5)              Core oversampling factor
        (mode = "ql")
GIMFIT2D is the heart of the GIM2D package. This is the task that actually performs bulge/disk decompositions. Therefore, it will require significantly more explanations than the others.

The Metropolis Algorithm

Historical trivia: Nicolas Metropolis was a mathematician who worked on the construction of the Los Alamos MANIAC computer used for hydrogen bomb calculations. Edward Teller, the so-called father of the H bomb, is one of the co-authors on the original paper (Journal of Chemical Physics, 1953) describing the algorithm. I bet you never expected to use a piece of the original H-bomb research in your own work. I sure did not.

The galaxy 2D light model used by GIMFIT2D currently has eleven parameters: total flux, bulge fraction (0=pure disk system), bulge effective radius re, bulge ellipticity e, bulge position angle, exponential disk scale length, disk inclination (0=face-on), disk position angle, subpixel x and y offsets of the galaxy center and finally, the background level. This 11-dimension parameter has a very complicated topology with lots of local minima in many cases. It is therefore important to choose a method which does not easily get fooled by those local minima. The Metropolis algorithm is such a method. It should be said right away that it is NOT an "efficient" algorithm i.e. it is very CPU intensive. A gradient search is efficient. On the other hand, a gradient search is "greedy". It will start from the initial parameter values, dive in the first minimum it encounters, and claim that the minimum is the global one. The Metropolis tries to do the opposite as explained below!

The Metropolis starts from an initial set of parameters given by the user and computes the likelihood P_0(w|D,M) that this set is the right one i.e. the one Mother Nature used. P(w|D,M) is the probability that the parameter set w is the right one given the data D and the model M. It then generates random perturbation about that initial location with a given "temperature". If the search is hot, large perturbations are tried. If the search is cold, small perturbations are tried. After each trial perturbation, the Metropolis computes P_1, the value the likelihood function at the new location, and immediately accepts the trial perturbation if P_1 > P_0. On the other hand, if P_1 < P_0, the Metropolis will accept the trial perturbation only P_0/P_1 of the time. This means that the Metropolis will sometime accept trial perturbations which take it to regions of lower likelihood. Why is this strange behavior so valuable? Well, if the Metropolis finds a minimum, it will try to get out of it, but it will have a finite probability (inversely proportional to the depth of the minimum) of succeeding.

In the current implementation of the Metropolis, convergence is achieved when the difference between two likelihood values separated by 100 iterations is less than 3-sigma of the likelihood value fluctuations.

Another important aspect of the Metropolis algorithm is its treatment of confidence intervals. In a typical fitting algorithm, formal fitting errors are derived from a curvature matrix computed around the best solution. This is perfectly reasonable when the maximum around the best solution is symmetric. This is not always the case especially at low S/N ratios where the maximum may be broad and asymmetric. In this case, the curvature matrix will not account for the presence of tails in the distribution. After the Metropolis converges, it Monte-Carlo samples the region where the likelihood is maximized and stores the accepted parameter sets as it goes along to build the distribution P(w|D,M). Once the region has been sufficiently sampled, the Metropolis computes the median of P(w|D,M) for each model parameter as well as 99% confidence limits. Note that this is different from least-squares fitting which would choose the maximum of P(w|D,M) as its best solution, a choice which is less robust to noise than the median value. This is shown in gory details in this section of my Ph.D. thesis on the suspicious coin experiment. In the limit of high S/N ratios, least-squares and the current Bayesian approach give the same answers.

Furthermore, the sampling of parameter space in GIMFIT2D shapes itself to the local topology. Every 50 accepted iterations, GIMFIT2D computes the covariance matrix of those iterations and inverts this matrix through Choleski inversion to produce a "step matrix". This step matrix is then used to generate step sizes for the next iterations. Under the control of the step and covariance matrices, the volume sampled by the trial steps continuously changes shape to respond to local conditions. I like to refer to that sampling volume as the "multidimoid". In plain english, if GIMFIT2D finds itself in "valley" of parameter space, it will start poking more and more along the axis of that valley in order to get out. There is no point for the multidimoid to try to climb steep walls when it might be easier (i.e. more probable) for it to simply walk out. For more information, see Vanderbilt and Louie, Journal of Computational Physics, 56, 259 (1984).

What does GIMFIT2D do? Description of the GIMFIT2D parameters

in_image is the name of the image for which a bulge/disk decomposition must be performed. This image should have odd dimensions. It is usually an postage-size image extracted from the main science image with XGAL.

out_image is an image of the PSF-convolved best model found by GIMFIT2D. GIMFIT2D automatically subtracts out_image from in_image to create the residual image res_image.

sym_image is the symmetrized version of in_image. Created only if dosym=yes.

msk_image is the postage-size section of the SExtractor segmentation image corresponding in_image. It is usually extracted by XGAL at the same time as in_image. The key to the flag codes recognized by B>GIMFIT2D are: The -1 flag is for cases where you did not have to run SExtractor to deblend galaxy images. So, when you do not use a segmentation image, msk_image should simply be the name of the image which have all its pixels set to -1.

sig_image is the image specifying the 1-sigma weight of each pixel for the fits. Set to "none" to use GIM2D's simple internal noise model.

logfile is the name of the output used by GIMFIT2D to store information about all the iterations tried by the Metropolis. Every line begins either one of the following letter codes: B, F, I, R, A, E, and P. The key to these codes is: There are 2 lines in logfile per iteration. The first line gives iteration number, number of iterations accepted so far, ln(likelihood), and parameter values tried in the current iteration. The second line gives the sizes of the 12 parameter perturbations tried in the current iteration.

Finally, at the very end of logfile, there are 16 lines giving information about the residual image indices, concentration and asymmetry indices. The CHI line gives the reduced CHI2 of the best model fit. See GREDCHI2 for more details. The RHALF line gives the half-light radius of the object computed by numerical integration of the best structural parameters. This radius is therefore the intrinsic radius of the object before PSF convolution. The RT, RA, AZ and DZ lines give estimates for 4 different asymmetry indices computed at various distances from the center of the objects. See the GRART page for details. The C1A-C4A lines give Abraham's concentration index (1st number, alpha = 0.1 to 0.4), Abraham's background-corrected asymmetry index (2nd number) and the background correction applied to the asymmetry index (3nd column).

mdparfile is a file which specifies the dimensions of the search volume to be explored by the ICF. Click here for a example of an .mdpar file. This file must contain 13 lines. The first line specifies the number of models the ICF has to create. The remaining twelve lines give set the initial conditions for the 12 galaxy model parameters. These lines are ordered as follows: L_tot, b, re, e, phi_b, r_d, i, phi_d, dx, dy, db and nb. Each line should have five numbers: the first two should be set equal to one another and should be equal to the initial parameter value. The third and fourth numbers should be the minimum and maximum hard limits on the parameter value. No values outside those hard limits will be explored by the ICF and the Metropolis. Finally, the fifth number is the initial size of the perturbation applied on the parameter. Say you want to explore bulge fraction between 0.4 and 0.8, the line in mdparfile for that parameter would read: 
0.6 0.6 0.4 0.8 0.4
I encouraged you not to be shy. Set really wide limits in mdparfile if you have no a priori knowledge of the structural parameters of your objects of interest. In fact, you have no choice since you better make sure your bulge fraction is allowed to freely go from 0 to 1 if you are trying to classify the morphologies of a sample of galaxies. On the other hand, let us say that you would like to "freeze" a parameter to its initial value. Why? Well, you might want to fit a pure bulge or a pure disk model to a particular galaxy. In the case of a pure disk fit, the line in mdparfile for the bulge fraction would read: 
0.0 0.0 0.0 0.0 0.0
In another example, say you would to freeze the disk position angle phi_d to +65 degrees. Then, the line in mdparfile for the disk position angle would read: 
65.0 65.0 65.0 65.0 0.0
If initparam is set to yes, then some initial parameter values from mdparfile will superseded by values derived from image moments. New initial values will be re-calculated for all parameters except the bulge fraction. The sky background level and sigma (given by bkgsig) will also be both superseded if dobkg is set to yes.

imscale is the image scale in arcseconds per pixel. Only useful if psftype is set to gaussian.

psftype is the type of Point-Spread-Function (PSF) to use. There are three options: gaussian, tinytim and user. If psftype is set to tinytim or user, then the PSF path and name should be given by in_psf. If psftype is set to gaussian, then make you provided a FWHM estimate with the parameter seeing. If you are using psftype=user, then your PSF image should have odd dimensions, and the PSF should be centered at ((nx-1)/2+1, (ny-1)/2+1) within that image

c_abs is the disk internal absorption coefficient. Values can range from 0 (optically thin) to 1 (optically thick).

bkgsig is the local background sigma either calculated from first principles if you know your detector and your reduction pipeline well, or directly derived from image statistics.

nsamp is an important control parameter for the Metropolis parameter. After achieving convergence, the Metropolis will have to accept nsamp iterations before it can stop. Thus, nsamp is the number of Monte-Carlo samples that will go into the calculation of the best parameter values and their respective confidence limits. There is no equation that will give you the right value for nsamp. You have to rely on intuition and experience. For example, if you are dealing with relatively high S/N ratio objects, the maximum likelihood region in parameter space should be relatively narrow and symmetric, so nsamp could be set to a relatively low value. On the other hand, if your objects have either low S/N ratios or are small, then the maximum should be shallow and asymmetric, so you will probably have to use a higher value of nsamp to achieve the same parameter space sampling density. I use a value of 300 for nsamp.

ccdgain is obviously the gain in electrons/DU of your detector. However, I have started using it in a slightly different way. I have taken the habit of keeping nframe set to one and defining ccdgain as the effective gain of the combined image stack. The effective gain will be different depending on whether the images have been averaged together or summed. If the detector gain is g and N frames were averaged together, then ccdgain should be set to N*g. If N frames were summed together, then ccdgain should be set to g. Always keep nframe equal to one for now. I will get rid of it in the near future.

metseed is a seed for the random generator producing Metropolis perturbations.

osamp is the oversampling factor for the core of the model. For the Wide Field camera, osamp should be 5. For the Planetary Camera, osamp should be 3. For ground-based data, osamp should be 1.

Last Update: 21 September 2000
Comments/questions to Luc.Simard(at)nrc-cnrc.gc.ca