COPYRIGHT (c) 1999
Kapteyn Astronomical Institute
University of Groningen, The Netherlands
All Rights Reserved.
Program: XGAUFIT
Purpose: XGAUFIT estimates parameters for multi component one
dimensional gauss functions using data in selected
profiles. It can use these estimates to fit the
parameters to the data in a least squares fit.
Output is send to subsets in a set.
Category: ANALYSIS, PROFILES, FITTING
File: xgaufit.src
Author: M.G.R. Vogelaar
(GUI: J.P. Terlouw)
Keywords: This is a list of keywords from an example run of the
application. It should be obvious to which gui elements
these keywords correspond.
AMPPERC= 50
AXUNITS= KM/S
BOXAXES= RA DEC
CRITAMP= 0
CRITDISP= 0
DATAUNITS= W.U.
FUNCTION= 0
GAUSET=fittedprofiles
GRIDRANGE= -30 30
LAB= 0.01
MAXITS= 200
NCOMP= 1
OUTSET=fittedparams
PARAM=
PHYSRANGE= 552.668 1548.98
PROFILE= freq
PROFINDX=
Q= 1 2 3 4 5
RMS= 1
SETNAME= real
TOLERANCE= 0.001
VSYS= 1050
Z0_EST= 0
Z0_FIX= YES
Z1_EST= 0
Z1_FIX= YES
Z2_EST= 0
Z2_FIX= YES
Buttons etc.:
Notes about the Graphical User Interface (GUI):
User input in this application is coupled to the Hermes keyword
mechanism. Keywords can still be specified on the Hermes command
line and also the syntax for input of numbers and text is the
same as for keywords. Also, Emacs-like text editing in input
fields is possible and we added mouse functionality.
E.g. the right mouse button will erase text in an input field.
With the middle button you can paste text stored in your mouse
buffer.
, and ctrl-U
Wipe current text from input field
ctrl-A Go to beginning of field
ctrl-E Go to end of field
ctrl-K Wipe from cursor position to end of
line (ctrl-Y does NOT yank back)
Go to next input field.
Control keys in the Help window:
ctrl-S Search forward
ctrl-R Search backwards
ctrl-V Scroll text one page down
(same as Page-Up)
ctrl-P Scroll text one line up
ctrl-Z Scroll text one line down
Notes:
1) Do not close pop-up windows using the cross symbol
on the upper windows bar. This will crash your application.
Instead use the close button or repeat pressing the button
in the main gui with which you opened your pop-up window.
2) On some systems, a menu item is un-selectable if the
Caps-Lock is on!
MAIN WINDOW:
Setname: input field.
Enter name of set which contains your profile data.
Profile axis: Menu
List of possible profile directions.
Grid range: Input field.
Default is the entire profile length in grids.
If you change one of the values, the corresponding
physical range is altered too.
Phys range: Input field.
Default is the entire profile length in physical coordinates.
If you change one of the values, the corresponding
grid range is altered too.
Box: Input field.
The limits in your set between which profiles are examined.
A box consist of two vectors. The first gives the low limits and
the second gives the high limits. The dimension of
a vector is the set dimension - 1. If you press
you will get a default box which includes all profiles in
the set.
Vsys: Input field
Value used in sorting when you want to sort in distance
to a position along the profile axis. Usually this will be
the systemic velocity of a galaxy.
Parameter output: Input field
Name of the output set on disk which holds all the fitted
parameters in subsets. Size of the output set is determined
by the program.
Fitted data output: Input field
Name of the output set on disk which holds the fitted data.
The output set is copied from the input set. Initially all
data is set to blank. You can get a residual cube by subtracting
this set from the input set e.g. with task SUB.
GO/STOP: Button
Start calculating initial estimates and fitting the
parameters to the data in the profiles. If you want to
abort this process, press this button again and wait a while
for the program to respond.
BUTTONS in button bar:
FILE
====
Menu options:
1) Get estimates from these subsets instead of calculating them.
You must give the name of the set and a subset specification.
E.g. myestset param 1 2 3.
You need as much subsets as components x 3.
2) exit program
FILTERS
=======
Amp. range: Input field
No values means no filtering. One value filters out all
amplitudes below this value and two values filter out all
amplitudes below the first and all amplitudes above the
second value.
Vel. range: Input field
Disp.range: Input field
See above. These last input fields need values in physical
coordinates.
FIT OPTIONS
===========
Tolerance: Input field
Stop criterion for least squares fit (see description below).
Lambda: Input field
Mixing parameter (see description below)
Max. its: Input field
If the stop criterion is not reached after this number
of least squares iterations, then abort the fitting of the
current profile. Write blanks in the output set in the
corresponding position.
RMS: Input field
Set the value with which the value of chi-squared is scaled to
get a value for the reduced chi-squared that is related to
the noise in your profile. The default is 1.
Min. ampl: Input field
Discard initial estimates with an amplitude lower than this
value.
Min. disp: Input field
Reject initial estimates with a dispersion lower than this
value.
Q: Input field
Smoothing factors (see description below)
PARAM
=====
Set fit parameters to fixed or free in the fit. If you want
to set a parameter to fixed, you have to enter a value
for this parameter.
For the background terms, an initial estimate must be
given in the input fields. The defaults for all the
background terms are zero.
FUNCTION
========
Menu with the following functions (see description below)
Gauss
Gauss Hermite (h3)
Gauss Hermite (h3,h4)
Voigt
NCOMP
=====
Menu with possible number of components per function.
SORT
====
Menu with sort options:
None
Peak,
Disp,
Vsys
Select sort option if you have selected more than one
component to fit and want to sort these components.
See also see description below.
HELP
====
Pop up this dc1 document
Description: BASIC INPUT
The decomposition of profiles into Gaussian components
can be a powerful method in the study of the kinematics
of neutral hydrogen. The name of the set for which you
want to study your profiles is entered in the main gui
of the program. The input field is called 'Setname'.
If a set is known, then the profile direction is selected
from a menu labeled 'Profile axis'. Usually this will
be the FREQ or VEL axis of your data set. Note that you can
select a profile direction as long as your input set has
more than one axis so that input is not restricted to
a three dimensional data cube. The length of your
profile axis must be greater than 1 and less than 4096.
The length of a profile is set in the input fields
'Grid range' and 'Phys. range'. A default is displayed
after the name of a set and the profile axis are known.
You can change these values either in grids or physical
coordinates. If you change one field you see that the
values in the other field will change also.
HINT: Inspect with XGAUPROF whether you can
exclude unwanted data by decreasing the grid range. However
if your grid range is too small you can run into trouble
while fitting base lines.
The 'Box' field is blank. Pressing results in
a default box which is the entire data range in grids.
Note that the box can be given in any coordinates as you
were used to in the non-gui GIPSY tasks.
In general, the program will work on profiles in any
direction and instead of velocity you could think of any
kind of physical coordinate. The program automatically
inserts the correct units in the header of the
output set.
PREPARING THE FIT
Estimates calculated by the program
-----------------------------------
You can fit the parameters of four different functions
to the profile data. The results are written to a
GIPSY output set on disk, each parameter in a different
subset per profile position. But the least squares fit
needs a good starting point for these parameters and
an important part of this program is dedicated to finding
good estimates by fitting second degree polynomials to the
profiles (See: Schwarz, 1968, Bull.Astr.Inst. Netherlands,
Volume 19, 405-413). An important role here play the so
called smoothing parameter 'Q'. The initial estimate
routine uses 2Q+1 points around a maximum to fit a second
order polynomial. The coefficient of the second-order term
of the polynomial is an approximation of the second derivative
of the observed profile. With a moments method
the center and dispersion are estimated. The maximum amplitude
is derived from the observed profile.
In practice there
is not one smoothing factor that gives the best (in terms
of chi-squared) estimates for all profiles. Therefore you
can enter in the 'FIT OPTIONS' menu a series of smoothing
factors ('Q'), which are all used in all profiles to get the
best estimates. You can exclude estimates by entering
values for a minimum amplitude ('Min.ampl') and the
minimum dispersion ('Min.disp'), but usually you leave
these fields blank and set a filter on the least squares fit
results.
Note that there is one special case build in the program.
If a profile contains blanks or zero-level values except
for one, then the estimate routine fails to return an
estimate. One can treat this peak as a peak with
known amplitude and center, and a full width at half maximum
of one pixel (which is automatically converted to a dispersion).
However, the least squares fit routine will not fit such
profiles because it converges to zero dispersion.
Estimates from an input set
---------------------------
Alternatively it is also possible to read estimates from
a set. This set could be an output set of a previous run of
XGAUFIT. If not, then you must be sure that the box of
your estimates set matches the box of your input set
and that the data has the same physical coordinates as
the set with the profiles.
Not only the set name must be entered, but also the subsets,
e.g.: myestimates param 1 2 3
With this option it is possible to separate the data with which
you create estimates, and create fits. For instance you can
create estimates with clipped or conditional transferred
data and use these estimates to fit parameters using the
original profile data in a second step. Such methods avoid
the noisy parts of your data and result in 'cleaner' output.
THE FIT PROCEDURE
-----------------
The actual fitting is a least-squares fit of a function
to a set of data points. The method used is described
in: Marquardt, J.Soc.Ind.Appl.Math. 11, 431 (1963).
It is a mixture of the steepest descent method and the
Taylor method.
Weights
-------
The fit routine uses equal weights for all data that
is not blank. Blank values get weight zero.
Tolerance
---------
The fit is controlled by three parameters which are set
in the 'FIT OPTIONS' window. Most important parameter is
the stopping criterion 'Tolerance'.
Fitting is stopped when successive iterations fail
to produce a decrement in reduced chi-squared less
than TOLERANCE= The value cannot be less than a
certain minimum as set by the system. This means that
maximum accuracy can be obtained with TOLERANCE=0.0
Usually only signal is present in a small region in the
profile while the chi square is calculated over the entire
baseline. In this case the Chi square primarily represents
the fitting of the baseline instead of the gauss and then
a very tiny decrease of the Chi square might very well be
significant for the fitting of the gauss to the signal.
So one should be very careful not setting TOLERANCE= too
high. If TOLERANCE=0 then the minimization is stopped at
machine precision. Increasing the tolerance results in
a shorter process time. If the tolerance cannot be reached
within the maximum number of iterations, the fit failed and
blanks are written in the output set.
HINT: In most cases the default tolerance 0.001 gives fast
performance and good fit results.
Mixing factor
------------
Lambda in the 'FIT OPTIONS' window is the mixing parameter in
the least squares fit routine.
Lambda determines the initial weight of steepest descent
method relative to the Taylor method. It should be
a small value (i.e. 0.001).
Maximum number of iterations
----------------------------
Fitting is also controlled by the maximum number of
iterations 'Max. Its.' in the 'FIT OPTIONS' window. If the
maximum number of iterations is exceeded, then the fit
failed and blanks are written in the output set.
Fixed or free parameters
------------------------
In the 'PARAM' window you can alter the status of the
fit parameters. Default, all function related parameters
are 'free' parameters (i.e. to be fitted) and the
background parameters (a constant, a linear and a quadratic
term) are fixed at a value 0.0.
Error on fitted parameters
--------------------------
The errors on the fitted parameters are calculated as
real errors by means of the covariance matrix in the fit
and the average deviation between data and fit over the
entire profile. However, if the region with signal is
short compared to the full length of the profile, the
average deviation between data and fit is primarily
determined along the part of the baseline where no signal
is present and then it becomes almost equal to the noise
in the profile and the real errors determined in this way
virtually become formal errors.
Note that if Hanning smoothing has been applied, the subsets
are two by two correlated. The errors on the parameters
are then a factor sqrt(2) higher.
MULTI COMPONENTS AND SORTING
----------------------------
If you fitted parameters of multicomponent functions,
then it could be useful to sort the fitted values. Usually
you will sort components in INCREASING order of distance
to a user given value on the profile axis (often the
systemic velocity) as given in Vsys in the main window.
The sort menu is found in the upper right area of the
main gui. Other sorting options are: sort components in
DECREASING order of peak value and sort components in
DECREASING order of dispersion
FILTERS
-------
The program has three filters for the output. You can
enter ranges for amplitude, velocity (or whatever is
on the profile axis) and dispersion. Fitted parameters
outside these ranges are filtered out, i.e. they are set to
blank in the output. If no values are entered then nothing
will be filtered. If only one value is entered for one
of the filters, then everything BELOW this value is
filtered out. In the gui log you can monitor how many profile
fits are rejected by these filters.
Warning: The velocity filter is straightforward for standard Gauss
functions. If one fits parameters of the Gauss-Hermite
series, then the filter applies to parameter b which is
only the mean (velocity) of the profile distribution if the
parameters h3 and h4 are zero.
FUNCTIONS
---------
The function menu provides the choice of four functions
of which you can fit the parameters to the data in your
profiles.
1) A standard Gaussian. Parameters are amplitude, center,
and dispersion,
2) Gauss-Hermite polynomial (h3). Parameters are a,b,c,
(which are NOT the same as amplitude, center and
dispersion) and h3. The parameters are translated to
the more familiar amplitude, center, dispersion and
skewness. All parameters mentioned above are saved
in the output set.
3) Gauss-Hermite polynomial (h3, h4). Same as above, but
an extra parameter h4 is included. The value for this
h4 parameter is translated to kurtosis.
Note that if one of h3 and h4 is not zero, the mean of
the distribution is not the position of the maximum.
The program calculates this parameter separately.
(Reference; Marel, P. van der, Franx, M., A new method
for the identification of non-gaussian line profiles
in elliptical galaxies. A.J., 407 525-539, 1993 April 20).
4) Voigt. Fit Lorentzian and Doppler widths together
with center and amplitude.
The mathematics behind all this is worked out
in a separate document:
(http://www.astro.rug.nl/~gipsy/xgauprof/index.html)
OUTPUT
------
A name is sufficient to specify the output set.
If the set exists, you are asked to confirm to write
it over. Before writing any data, an existing set will
always be deleted first to avoid complaints about
incompatible sizes if you change boxes or grid ranges and
do not alter the name of the output set.
The program creates the new set on disk and copies relevant
header items. It replaces the profile axis
by a new axis called PARAM-XGAUFIT. Its length
depends on the maximum number of components to be fitted and
the selected function. It is also possible to save the
estimates only (Est.only button in 'FIT OPTIONS' window).
The Hermes log displays the structure of your output.
This information, together with relevant units, is also
written to the comment fields header of the output
set so that you can always retrieve subset information
e.g.:
HISTORY ITEM=C MODE=R INSET=myout
or:
FIXHED INSET=myout ITEM=
which shows header items at set level or
FIXHED INSET=myout p ITEM=
to show units , minimum, maximum and number of blanks at
subset level.
Not only the fitted parameters are saved,
but also the errors on these parameters, the background
parameters, total flux, number of iterations, reduced
Chi-squared of the fit, the standard deviation of the residuals
and the smoothing factor 'Q' used by the estimate routine.
After finishing the fitting process, the program writes
a description of the subsets to the comment fields of the
header of the output set. The units, minimum and maximum values
and the number of blanks are written on subset level in
the FITS keywords BUNIT, DATAMAX, DATAMIN and NBLANK.
Output of the fitted profiles
-----------------------------
There is an option to reconstruct your data using the fitted
parameters. These 'fitted profiles' are written to a set
which is a copy of the input set, i.e. the header is a copy
and the data is set to 'blank'. The size of the output is
limited by the values for the box and the profile range
(Grid range or Phys. range). The header keywords for the
minimum and maximum data value and the number of blanks are
NOT updated. This update can be done with task MNMX.
Each fitted profile is written to this output set. If the fit
was not successful a blank appears in all the profile subsets
of the output set.
A message is appended to the HISTORY of the output set.
The message contains the task name and time/date.
Updates: Jul 20, 1999: VOG, Document created.
Jul 27, 2000: VOG, g->Yfit was not initialized to NULL.
This seemed only to cause problems on PC's(?)
Nov 6, 2001: JPT, Wait for processing finished before quitting.
Apr 15, 2009: VOG, Removed macro NINT
Jun 28, 2010: JPT, Changed log font.