Table of Contents
List of Equations
Table of Contents
Fityk is a flexible and portable program for nonlinear fitting of analytical functions (especially peak-shaped) to data (usually experimental data). In other words, for nonlinear peak separation and analysis.
It was developed for analyzing diffraction patterns, but can be also used in other fields, since concepts and operations specific for crystallography are separated from the rest of the program.
Fityk offers various nonlinear fitting methods, background subtracting, various data transformations, easy placement of peaks and changing peak parameters, automation of common tasks with scripts, and much more. The main advantage of the program is a flexibility - parameters of peaks can be arbitrarily binded with each other, eg. width of peak can be an independent variable, can be the same as width of other peak or can be given by complicated - common for all peaks - formula.
The GUI was written with portable wxWidgets library, so fityk works in both Unix/Linux and MS Windows environments.
As far as I know, the only program which also is specialized in fitting peaks is PeakFit. Some of useful features found in PeakFit, that are not present in fityk yet, will be implemented in future.
Fityk is free software; you can redistribute and modify it under the terms of GPL, version 2. See Appendix E, License for details. You can download the latest version of fityk from http://www.unipress.waw.pl/~wojdyr/fityk . or http://fityk.sf.net . To contact author, visit the same page.
After this introduction, read entire Chapter 2, Using fityk. It is also recommended to read Chapter 3, Command reference , especially the section called “General syntax ”, the section called “Sum of fitting functions ” and, if your field is crystallography, reading the section called “Module for crystallography” is a must. Note, that menu commands and toolbar buttons are generally not explained, because they are either intuitive and self-explaining or analogical to textual commands, that have accurate description.
Table of Contents
The program comes in two versions: GUI (Graphical User Interface) version - more comfortable in most cases, and CLI (Command Line Interface) version (named cfityk to differentiate, Unix only) - used only in special situations, when visualization is not essential and user prefers to use a text terminal.
If CLI version is compiled with the GNU Readline Library, it enables command line editing and command history like in bash. Especially useful is TAB-expanding. All fityk commands, filenames and option names can be expanded. Options are set by d.set, s.set, f.set etc. commands. In CLI version, data and sum of functions can be visualized with gnuplot program. It is done automatically, if gnuplot is found.
GUI version is written using wxWidgets library. One of the main features of this library is portability. Program can be run on most Unix species with GTK+ (it is developed on Linux) and on MS Windows. Ports to other platforms supported by wxWidgets library are also possible. Read the section called “Graphical interface ” to learn how to use it.
Let us analyze a diffraction pattern of NaCl. Our goal is to determine the position of the center of the biggest peak. It is needed for calculating pressure, under which the sample was measured, but the further processing does not matter in our example.
Data file used in this example is distributed with the program and can be found in samples directory.
First load data from file nacl01.dat. You can do it by typing d.load 'nacl01.dat' in CLI version (or in GUI version in the input box - at the bottom, just above the status bar). In GUI, you can select -> from menu and choose proper file, instead. All commands can be abbreviated with letter-dot-letter, eg. you can type: d.l 'nacl01.dat'.
If you use GUI, you can zoom in the biggest peak using Shift and button on the auxiliary plot. Other ways of zooming are described in the section called “Mouse, Modes and Toolbars ”. If you want the data to be drawn using bigger points or line, or if you want to change color of line or background, press mouse button on main plot and use or menu from pop-up menu.
mouse button on the auxiliary plot (the plot below main plot). To zoom out, press either button orNow all data points are active. Because only the biggest peak is our interest, the rest of points can be deactivated. Type: d.range 23.0:25.0 or change to range mode (press Data-Range Mode button on toolbar) and select range to be deactivated with mouse button.
We see that our data has no background, we would have to worry about, so now we only have to define peak with reasonable initial values of peak's parameter and fit it to data. We will use Gaussian. To see it's formula, type: s.info ^G or look for it in the documentation (in Appendix C, List of functions).
To define peak, type: s.add ^G ~60000 ~24.6 ~0.2 or select Gaussian from the list of functions on toolbar and press auto-add toolbar button. There are also other ways to add peak in GUI, try add-peak mode.
Now let us fit the function. Type: f.run or select
-> from menu.To see position of center, type: s.info ^0 ('^' stands for function and ^0 for first function) or s.info @1 ('@' stands for simple-parameter, @0 is the height and @1 is the center). This will be explained later in this manual. In GUI information about peak is displayed on status bar, when you move mouse over top of peak.
That is all. To do the same second time, you can write all the commands to file, eg. to file nacl01.fit, and use this script: o.include 'nacl01.fit' or select
-> from menu, or run program with the name of script: bash$ fityk nacl01.fitThe window of fityk program consists of (from the top): menu bar, toolbar, main plot, auxiliary plot, output window, input field and status bar. The input field allows to type and execute commands in similar way as it is done in CLI version. The output windows shows results. Click the mouse button on it to see its configuration options.
Most of menu commands correspond to commands described in Chapter 3, Command reference . Some of them are simply translated to text strings and executed - you can see it in output window. I hope that menu commands are self-explaining (especially after reading the rest of this manual) - just try them.
The main plot can display data points, functions, phases and/or the sum of all functions. Use pop-up menu (click Phase is a sum of related functions. The word comes from the crystallographic module, where a phase of investigated powder is represented by a set of peaks.
button on the plot) to configure it.The auxiliary plot can display a difference of data and the sum of functions, the difference divided by standard deviation of data or another plot. It can be controlled by its pop-up menu. I hope having a look at this menu and a minute of experiments will explain possibilities of the auxiliary plot.
Sometimes it is useful to print effects of our fitting work. Hard copy will contain plots visible on the screen, scaled to fit the page. The only difference is that backgrounds of plots will be white (to not waste toner/ink). So it may be necessary to change color of axis or data to darker one.
Configuration of GUI (everything that can be changed using pop-up menus and size of the window) can be saved using
-> . Two different configurations can be saved, what allows easy changing of colors for printing. On Unix platform, these configurations are stored in file in user's home directory. On Windows - they are stored in registry (perhaps in future they will be also stored in file).Now it is a time to explain quite complicated mouse usage. The usage of the mouse on menu, dialog windows, input field and output window is intuitive, so the only topic described here will be: how to effectively operate mouse on plots.
Let us start with the auxiliary plot. Shift) will zoom it out and display all data.
button displays pop-up menu, with button you can select range to be displayed - range on x axis. Clicking with button (or with button and pressedOn the main plot, meaning of the Ctrl (or Alt or Shift). button can be used to select a rectangle, that you want to zoom in. If an operation has two steps, like rectangle zooming (first you press button to select first corner, then you move mouse and release button to select second corner of rectangle), you can cancel it by pressing another button when first one is pressed.
and mouse button depends on current mode, that can be changed using toolbar or menu. You will see hints on status bar. In normal mode, button is used for zooming and invokes pop-up menu. The same behaviour can be obtained in any mode by pressingTable of Contents
All commands begin with a token x.xxxxx (letter-dot-letter). Character before dot denotes what the command is related to, and the word after dot describes what it does. Every command can be abbreviated with x.x (letter-dot-letter). The letter before dot can be one of:
data related commands
sum related commands
manipulations, various algorithms and computations related commands
fitting related commands
other commands - including scripts, logging to file, plotting etc.
commands from crystallography/diffraction module
For instance, f.run is used to start fitting and f.info to display some informations related to fitting.
Commands are finished with new line (or by pressing enter interactively), semicolon (;) or comma (,). Basically, one command is in one line. If, for some reasons, it is more comfortable to place more than one command in one line, they can be separated with semicolon (;). Sequence of commands with the same beginning token can be replaced with one command, which has one common beginning token and tokens of each command separated by commas, eg.
s.info ^0, ^1, ^2
Symbol '#' starts a comment - everything from hash (#) to the end of the line is ignored.
There is one command, that has the same syntax in all families of commands. It is x.set command, where x is one of letters described above. It is used to set or read values of various options. Command
shows the list of all options and their values if option is not given, or only value of option. Command:
is used to set a new value of option. In CLI version, TAB can be used to expand names of options and possible values. In GUI version, dialogs can be used.
Data are stored in files. Unfortunately, there are various formats of files with data. Usually we have text files and one line in file corresponds to one data points. Every line should contain at least two numbers: x and y of point. It can also contain standard deviation of y coordinate. Numbers can be separated by whitespace, commas, colons or semicolons. Some lines can contain comments or extra informations. If these lines have a hash (#) in first column, they are ignored. In other case, they are usually also ignored. There are also other file types, that can be read: .rit and .mca. In future, the way the special file formats are handled will be changed.
Points are loaded from files using command
where xcol, ycol, scol, are unsigned integers. Note that the name of the file is inside of single quotation marks. If the file is in a normal format (columns with numbers) it can be specified which column contains x, y and, optionally, std. dev. of y.
Command
can export data to ASCII file, that has 3 columns: x, y and standard deviation. Only active points are being exported (see next chapter to learn about active and inactive points). + will cause, if specified file exist, appending to the file rather than overwritting.
Some information about current data can be obtained using command:
We often have situation, that only part of data from file is interesting for us. We should be able to exclude selected points from fitting and all computations. It can be done with command d.transform (which will be explained later) or with mouse-click in GUI. The idea of active and inactive points is simple: only the active ones are subject to fitting and peak-finding, inactive ones are neglected in these cases.
When fitting data, we assume that only y coordinate of data is subject to
statistical error in measurement. Is is very common assumption.
To see how y standard deviation
influences fitting (optimization), look at
weighted sum of squared residuals formula in the section called “Nonlinear optimization ”.
We can also think about weights of points -
every point has a weight assigned, that is equal
Standard deviation of points can be read from file together with x and y coordinates. Otherwise, it is set to max(sqrt(y), 1.0), Setting std. dev. as a square root of value is common and has theoretical grounds when y is the number of independent events.
In fityk 0.5.0 new powerful command d.transform was introduced. It's not documented here yet. The best documentation for a while is the extensive comment at the beginning of datatrans.cpp source file. You can see CVS version of the file here: http://cvs.sourceforge.net/viewcvs.py/fityk/fityk/src/datatrans.cpp?only_with_tag=HEAD&view=markup
Historically, there were separete commands for dealing with spline background and with calibration. Now all changes to data points are made using d.transform command, but there is still background mode in GUI, which helps in manual baseline substraction and generates proper d.transform command to apply the change to the data. Have a look at two useful GUI options related to background mode: background type and minimal BG points distance.
Sometimes the instrument introduces errors to the x coordinate of the data points, which can be determined by using a reference sample. Correcting these errors is called here calibration. (See also description of zero-shift.) Historically, there was a special command for calibration in fityk, but it was removed to simplify the program. Now you can change x coordinate of data points using d.transform command, analogically to baseline substraction.
Let call a set of data that usually comes from one file - a dataset. All operations described above concern one dataset, but you can also work with many datasets. Datasets can be grouped in plots. There can be any number of plots, and any number of datasets in each plot. There is always one active plot and one active dataset in each plot. Inactive datasets are drawn, but do not influence any operations. So there are only two reasons to have more than one datasets in one plot - it is possible to compare them visually and to easily switch between them. Let me repeat: all command are related only to the active dataset.
A so-called plot contains not only datasets, but also described below so-called sum. Only datasets and sum in active plot are drawn and all commands except fitting are related to active plot. But for f.xxx commands it doesn't matter if plot is active - all plots are fitted simultanously. This is why one may want to have many plots - to fit simultanously datasets with functions, that can share its parameters.
You can manage plots and datasets using the data pane in GUI or d.activate command. The following syntax:
is used to switch to plot p and dataset d. It is possible to change only active plot or to change only active dataset. To create (and activate) new plot or new dataset replace p or d with asterisk (*).
The sum of functions S - curve that is fitted to data - is itself a
function. The value of the whole sum is computed as a sum of functions,
like Gaussians or polynomials, and can be given by formula:
,
where fi is a function of x, and
depends on a vector of parameters a. This vector contains all
fitted parameters.
Because we often have the situation, that the error
in the x coordinate of data points can be modeled with function z(x; a),
we introduce this term to sum:
where
. Note that the same z(x) is used in all functions.
Now we will have a closer look at fi functions. Every function fi has a type chosen from function types available in program. One of these types is the Gaussian. It has the following formula:
a0 is a height of peak, a1 - position of center and a2 - half width at half maximum. Every of these parameters can be:
a constant parameter - a constant real number.
a simple parameter - a real number, that can be changed.
compound parameter - a function of parameters, eg. sum of two simple parameters or product of constant and simple parameter. I called it g-parameter (this name will be changed to more meaningful in future). A g-parameter is of course different from functions fi. The basic difference is that g-parameter is not a function of x.
Functions fi are called f-functions. See the list of available f-function types in Appendix C, List of functions.
Functions zi(x;a) are called zero-shifts. This name is not precise (FIXME: what is a better name? calibration function? instrumental function?). The simplest of them has a formula: z = a0 and can be used for modeling error of shifting zero in instrument. In general case, so-called zero-shifts are parametrized functions of x, like f-functions.
Every f-function, zero-shift and g-parameter has its type and its parameters - the number of its parameters is determined by the type. Two examples:
Gaussian is a type of f-function, and f-function of this type has 3 parameters: a0, a1 and a2 (see formula above).
One of g-parameter types is division: a0/a1. A g-parameter of this type contains two parameters.
F-functions, zero-shifts, g-parameters and simple-parameters are elements of a sum. Every sort of element has a symbol: f-function (^), zero-shift (<), g-parameter ($) and simple-parameter (@). Every f-function contained in the sum has a number assigned, that allows referring to this f-function. If there is only one f-function defined, it has number 0. If there are two f-functions, the first has number 0 and the second 1. As you will see further in this chapter, commands refer to f-functions using symbol (^) and number, eg. ^0 or ^15. The same about other sorts of elements: first zero-shift is referred to as <0, fifth simple-parameter - as @4, etc.
Three kinds of parameters were described above. They are used in f-functions, zero-shifts and g-parameters, and are denoted in adequate formulae as ai. Note, that actually only simple-parameters are fitted. The term parameter used in context of fitting and the term vector of parameters refer to simple-parameters only.
Every simple parameter has a value and, optionally, domain. Command
creates simple-parameter equal p. Once created, the parameter will be given a number (first available), say m, and later can be referred to as @m. Domain is used only by fitting algorithms, when they need to randomly initialize variable, or change it. In some cases only width of domain is needed. Domain must be enclosed by square brackets. It can be specified in three ways:
[a:b] -- from a to b,
[a +- b] -- from a-b to a+b,
and [ +- b] -- from v-b to v+b, where v is the value of simple parameter.
If domain of parameter is not defined, it is assumed to be [v - r v, v + r v], where v is the value of simple parameter and r is the value of option default-relative-domain-width.
F-functions, zero-shifts and g-parameters can be created using command:
where t is a letter denoting type of f-function, zero-shift or g-parameter (respectively ^t, <t or $t). It is followed by parameters in number required for given type of object. Parameters can be given as:
@n - existing simple parameter, eg. @15,
$n - existing g-parameter, eg. $15,
_p - constant value, eg. _3.14159,
~p - simple parameter created on-the-fly, eg. ~3.2
In particular, a g-parameter can be defined as a function of another g-parameters. As a result of it, every parameter of f-function (eg. of Lorentzian function) can be an arbitrary combination of simple parameters. You can learn about available types using command s.info described below.
Every f-function, zero-shift, g-parameter or simple-parameter, that was created (added), can also be removed. The only condition is, that there are no references to it. It is obvious - you can not remove a simple parameter, that represents a height of a peak, before removing the peak. The syntax of command is straightforward:
n is the number of f-function, zero-shift, g-parameter or simple-parameter. It can be replaced with an asterisk(*). Command s.remove ^3 will remove ^3, and s.remove ^* will remove all f-functions. To avoid gaps in numeration, all larger numbers are decreased by one. If f-function ^3 is removed, ^9 will become ^8 and ^4 will become ^3. All references are properly updated. If option recursive-remove is set, also parameters, in our example parameters of ^3, will be removed.
Information about every element of the sum can be obtained using command:
And informations about available types of f-functions, zero-shifts and g-parameters can be displayed with:
t is a letter that stands for type, eg. ^G means f-function Gaussian. If type is not specified, list of all available types will be displayed.
The whole sum or selected f-functions can be saved to file in various formats with command:
where t is a letter that stands for format of file: f - mathematical formula, d - xy points, p - peak listing, x - XFIT peak listing (used by some crystallographic software), s - fityk script. + will cause, if specified file exist, appending to the file rather than overwritting.
Value and domain of simple parameters can be changed with command:
p is a new value of @n, or, if % is given, new value is equal to p% of current value. domain is specified in the same way as in command s.add. If domain is not given, it is not changed.
Vector of simple parameters is a subject to change in a fitting process. Sometimes it is useful to refine only selected parameters. The parameters that are not to be refined are called frozen. You can freeze and unfreeze parameters with command:
This command will freeze @n if exclamation mark (!) is not present and unfreeze if it is present. To freeze or unfreeze all parameters, use asterisk(*) instead of n. Example: s.freeze @*, ! @0, ! @3 causes, that only @0 and @3 would be fitted. To display information about currently frozen parameters, just type:
Every change of values of parameters is remembered. It is remembered since last adding or removing simple parameter operation. You can display history of changes with command:
and move to selected position in history with command:
Without sign - you will move to n-th position, with plus(+) or minus(-) - you will go n steps forward or backward. The history is similar to a history in web browser - if you move backward and then add a new entry, eg. by fitting, last entries of history will be erased. In some cases it can be reasonable to fit curve with one method, move back to initial parameters, fit with another method and compare results. To make this comparison easier, you can save result obtained with first method using command:
The command toggles saved and not saved state of history item, what means that using it second time on the same item will allow to delete it. If n is not specified, current history item is used.
Now it is time to explain how value of the sum and its derivatives are computed. You may skip it and read next chapter, if you are not interested.
All fitting methods involve computing value of the sum S(x) for x = x1 x2 ..., xN and such series of computations are repeated multiple times, for different parameter vectors. x1 x2 ..., xN are x coordinates of active data points. In typical cases N >> 1, therefore, when thinking about time consumed by computations, only operations executed for every point matter.
At the beginning of every series of S(x) computations values of g-parameters and some expressions in f-functions, that do not depend on x, are computed. Then S(x) value of every f-function and zero-shift is computed in every point. To speedup computations, tails of peaks, that are very small, can be neglected. This can be controlled with two options: cut-tails and cut-tails-level. If option cut-tails is set, at the beginning of each series of computations, for every peak, points, where the peak has value approximately equal to value of cut-tails-level option, are computed. Outside of these points, value of the peak function is assumed to be zero.
Most commonly used nonlinear fitting methods (in our program Levenberg-Marquard) are using first derivatives of functions. The program knows derivative formulae for all types of f-functions, zero-shifts and g-functions, but it can also use numerical differentiation. If is useful only for testing the program. Numerical differentiation is used when option numerical-d is set. If option numerical-d-both-sides is set, it uses formula df(x)/dx=(f(x+h)-f(x-h))/2h, otherwise df(x)/dx=(f(x+h)-f(x))/h. h is computed as value of option numerical-d-rel-change multiplied by x.
A value either of the whole sum or of a an element of the sum, and all appropriate derivatives can be viewed using command:
eg. s.value ^0 41.5. With asterisk(*), derivatives, computed symbolically and numerically are shown and compared, what is useful for program developers.
This is the core. We have a set of observations (data points), and we want to fit a model (sum of functions), that depends on adjustable parameters, to observations. Let me quote Numerical Recipes, chapter 15.0, page 656 (if you do not know the book, visit http://www.nr.com ):
The basic approach in all cases is usually the same: You choose or design a figure-of-merit function (merit function, for short) that measures the agreement between the data and the model with a particular choice of parameters. The merit function is conventionally arranged so that small values represent close agreement. The parameters of the model are then adjusted to achieve a minimum in the merit function, yielding best-fit parameters. The adjustment process is thus a problem in minimization in many dimensions. [...] however, there exist special, more efficient, methods that are specific to modeling, and we will discuss these in this chapter. There are important issues that go beyond the mere finding of best-fit parameters. Data are generally not exact. They are subject to measurement errors (called noise in the context of signal-processing). Thus, typical data never exactly fit the model that is being used, even when that model is correct. We need the means to assess whether or not the model is appropriate, that is, we need to test the goodness-of-fit against some useful statistical standard. We usually also need to know the accuracy with which parameters are determined by the data set. In other words, we need to know the likely errors of the best-fit parameters. Finally, it is not uncommon in fitting data to discover that the merit function is not unimodal, with a single minimum. In some cases, we may be interested in global rather than local questions. Not, "how good is this fit?" but rather, "how sure am I that there is not a very much better fit in some corner of parameter space?"
Our function of merit is WSSR - weighted sum of squared residuals, called also chi-square:
Weights can be are based on standard deviations,
.
You can learn why squares of residuals are
minimized eg. from chapter 15.1 of
Numerical Recipes. So we
are looking for a global minimum of chi2.
This large field of numerical research - looking for minimum or maximum -
is usually called optimization; it is non-linear and global optimization.
Fityk implements
three very different optimization methods. All are well-known and
described in many book.
You can switch between fitting methods using command:
where f is one of following lower case letters: g - Genetic Algorithm, m - Levenberg-Marquardt gradient method, s - simplex method. If f is not specified, current fitting method is displayed.
f.method affects the rest of f.xxx commands.
To run fit, use commands
and
The first initializes and runs fitting algorithm, and the second starts fitting without the initialization. n is the maximum number of iterations. All non-linear fitting methods are iterative, and there are two common stopping criteria. The first is the number of iterations. If it is not specified with command, value of option default-max-iterations is used. And the second is the number of evaluations of objective function (WSSR), specified by value of option max-wssr-evaluations. It is approximately proportional to time of computations, because most of time in fitting process is taken by evaluating WSSR. There are also other criteria, different for each method. Note that values of all options, even those common for all methods, are set separately for every method.
If you give too small n to f.run command, and fit is stopped because of it, not because of convergence, it makes sense to use f.continue command to process further iterations. [TODO: how to stop fit interactively]
When fit is running, each iteration outputs some informations. If they are scrolling too fast, you can reduce it n times by assigning n to output-one-of option.
Setting o.set autoplot = on-fit-iteration will draw a plot after every iteration, to visualize progress. (see autoplot)
Command [1] :
can be used to display WSSR, symmetric errors and variance-covariance matrix.
Available methods can be mixed together, eg. it is sensible to obtain initial parameter estimates using simplex method, and than fit it using Levenberg-Marquard method. Command s.history can be useful for trying various methods with different options and/or initial parameters and choosing the best solution.
Some fitting methods are using random number generator. In some situations one may want to have repeatable and predictable results of fitting, eg. to make a presentation. Seed for a new sequence of pseudo-random numbers is set at the beginning of fitting initialization (when f.run is called). If value of pseudo-random-seed option is set to -1, the seed is based on system time and sequence of pseudo-random numbers is every time different. In other case, if pseudo-random-seed option has a non-negative value, this value is used as a seed. Remember, that like all other options, value of pseudo-random-seed is independent for each fitting method.
It is a standard of nonlinear least-squares routines. It involves computing first derivatives of functions. For description of L-M method see Numerical Recipes, chapter 15.5 or Siegmund Brandt Data Analysis, chapter 10.15 (there is a Polish translation of the second). In a few words: it combines inverse-Hessian method (called Gauss-Newton method?) with steepest descent method by introducing lambda factor. When lambda is equal 0, the method is equivalent to inverse-Hessian method. When lambda increases, the shift vector is rotated toward the direction of steepest descent and the length of the shift vector decreases. (The shift vector is a vector that is added to the parameter vector.) If the better fit is found in iteration, lambda is decreased - it is divided by value of lambda-down-factor option (default: 10). Otherwise, lambda is multiplied by value of lambda-up-factor (default: 10). You can also change a value of lambda-starting-value option.
Marquardt method has one stopping criterium apart from common stopping criteria. If it happens two times in sequence, that relative progress in value of objective function (WSSR) is smaller then value of stop-rel-change option, fit is considered to be converged and is stopped.
L-M method finds a minimum quickly. The question is, if it is the global minimum. It can be a good idea to add a small random vector to the vector of parameters and try again. This small shift vector is added, when value of shake-before option is positive (by default it is 0). Value of every parameter's shift is independent and randomly drawn from distribution of type specified by value of shake-type option (see option distrib-type) in simplex method). The expected value of parameter shift is directly proportional to both value of shake-before option and width of parameter's domain.
This time I am quoting chapter 4.8.3, p. 86 of Peter Gans Data Fitting in the Chemical Sciences by the Method of Least Squares
A simplex is a geometrical entity that has n+1 vertices corresponding to variations in n parameters. For two parameters the simplex is a triangle, for three parameters the simplex is a tetrahedron and so forth. The value of the objective function is calculated at each of the vertices. An iteration consists of the following process. Locate the vertex with the highest value of the objective function and replace this vertex by one lying on the line between it and the centroid of the other vertices. Four possible replacements can be considered, which I call contraction, short reflection, reflection and expansion.[...]
It starts with an arbitrary simplex. Neither the shape nor position of this are critically important, except insofar as it may determine which one of a set of multiple minima will be reached. The simplex than expands and contracts as required in order to locate a valley if one exists. Then the size and shape of the simplex is adjusted so that progress may be made towards the minimum. Note particularly that if a pair of parameters are highly correlated, both will be simultaneously adjusted in about the correct proportion, as the shape of the simplex is adapted to the local contours.[...]
Unfortunately it does not provide estimates of the parameter errors, etc. It is therefore to be recommended as a method for obtaining initial parameter estimates that can be used in the standard least squares method.
This method is also described in previously mentioned Numerical Recipes (chapter 10.4) and Data Analysis (chapter 10.8).
After changing current fitting method to Nelder-Mead simplex, what can be done, as you already know, either using menu or using command f.method s you have a few options specific to this method, that can be changed. One of those is a stopping criterium min-fract-range. If value of expression 2(M-m)/(M+m), where M and m are values of worst and best vertices respectively (values of objective functions of vertices, to be precise), is smaller then value of min-fract-range option, fitting is stopped. In other words, it is stopped if all vertices are at almost the same level.
The rest of options is related to initialization of simplex. Before starting iterations, we have to chose set of points in space of parameters, called vertices. Unless option move-all is set, one of these points will be the current point - values that parameters have at this moment. All but this one are drawn as follows: each parameter of each vertex is drawn separately. It is drawn from distribution, that has center in center of domain of parameter, and width proportional to both width of domain and value of move-multiplier parameter. Distribution type can be set using option distrib-type as one of: uniform, Gaussian, Lorentzian and bound. The last one causes value of parameter to be either greatest or smallest value in domain of parameter - one of two bounds of domain (assuming that move-multiplier is equal 1).
Note: The title of this chapter is not appropriate now, but in future more commands will be implemented and described here.
For now, fityk offers only primitive algorithm for peak-detection. It looks for highest point in given range, and than tries to find out width of peak. This method is used by crystallographic module, and it is used when adding peaks (semi-)automatically in graphical version of the program. It can be called directly, using command:
If center and margin are not specified, highest peaks are searched in all data, otherwise only a range (center - margin ; center + margin) is considered. Default value of margin is equal to value of option search-width. Option estimate-consider-sum is set by default and causes that program is looking for peaks by analyzing difference of the data and the sum of functions. If not set, only data will be considered. Setting option cancel-peak-out-of-search disables adding peak, when the highest point in the range is found at the border of requested range.
Commands described in this chapter can not be classified into categories mentioned in previous chapters. Some of them are related to script files.
Command
reads and executes commands from file filename. If exclamation mark (!) is present, program is reseted to initial state. Only in this case filename can be omitted. You can specify lines or ranges of lines in file, that are to be read.
Command:
starts logging to file user input (i), output of program (o) or both of them (a). To stop logging, or to display informations about log file, use command:
With exclamation mark (!) - it will stop logging, otherwise - display informations.
There is also a command
that is writing a script, which can be used to restore current state of the program. It is different from logging, and usually output given by logging is easier to read by humans.
Command o.plot controls visualization of data and the sum. It is used to plot given area - in GUI it is plotted in program's main window, in CLI popular program gnuplot is used, if available. Actually there is no reason to type this command in GUI version, because it is very easy to zoom in and zoom out using mouse.
xrange and yrange have one of two following syntaxes:
{[} [min] : [max] {]}
.
The second is just a dot (.), and it remains appropriate range not changed. It will be explained on examples:
o.plot [20.4:50] [10:20] - plot will show x from 20.4 to 50 and y from 10 to 20 |
o.plot [20.4:] - x from 20.4 to end of data points, y range will be fitted to contain all data. |
o.plot . [:10] - x range will not be changed, lower bound of y range will be a bit smaller than the smallest y coordinate of point in this x range, and upper bound of y range is equal 10. |
o.plot [:] [:] - all data will be showed. |
o.plot . . - the visible area will not be changed, but information about this area will be printed. |
o.plot - the plot will be refreshed. |
There is also command in o.plot family, that decides about plotting background:
Plus (+) or minus (-) are causing background to be, respectively, added or not to the sum and the data while plotting.
Value of option autoplot changes automatic plotting behaviour. By default, plot is refreshed automatically after changing the data or the sum of functions. It is also possible to visualize every iteration of fitting method by replotting peaks after every iteration.
Value of verbosity option changes amount of text output. Possible values of verbosity and autoplot options can be expanded with TAB in CLI or are visible as a list in GUI.
Command:
makes program waiting sec seconds, doing nothing.
At last, if option exit-on-warning is set, any warning will close the program. It ensures, that no warnings can be overlooked.
This module simplifies using fityk for analyzing diffractograms. User must know approximate unit-cell parameters, and know or guess Miller's indices of reflections. Widths and shapes of peaks can be independend or angle-depended parameters. The whole pattern is refined. Note, that individual profiles can also be fitted, but this module is not needed for it.
It is assumed, that x coordinate of data points represent scattering angle
and wavelengths are constant.
Everything, what can be done with crystallographic commands, could be done with s.add, s.remove, s.info and m.findpeak commands, although it would be very inconvenient.
The wavelength can be specified with command:
lambda and intens are parameters and can be specified as parameters in s.add command. It practice, two forms are used: _p (eg. _1.93597) if wavelength is known or ~p (eg. ~0.9) if it should be fitted. intens is a relative intensity of wavelength, and it matters, when there are two or more wavelengths. Above command adds a wavelength with its ratio to a set of wavelength. This means, that if you want to define two wavelengths, you should use this command twice, and if you want to change constant wavelength, you must delete previously defined one and add a new one. To delete all defined wavelengths, use command
and to display informations about them, just type:
Now about a few essential options. You must choose a peak-shape, that will model a reflections. To do this, assign proper value to option peak-type. You can use Gaussian, Lorentzian, Pearson VII, pseudo-Voigt or Mod-TCHpV. The last one is also a pseudo-Voigt, but its width and shape parameters are given by following formulae:
Positions (centers) of peaks are determined by Bragg's Law, and heights (and, therefore, areas) are fitted almost freely - the only constraint is, that ratio of heights of peaks in one plane is equal to appropriate ratio of wavelength intensities. If you do not use Mod-TCHpV type, width and shape parameters of peak can be:
fitted autonomously for every peak,
equal for peaks in one plane (if there is more than one wavelength used, plane is represented by more than one peak),
equal for all peaks in one phase,
or function of
; for width, it is Lowe-Ma formula:
HK usually means FWHM, but note that above formula can be also applied to HWHM, by rescaling refinable variables: U'=U/4, V'=V/4, etc. Shape parameter of peak can be modeled by function of peak's position, for pseudo-Voigt peak:
and for Pearson peak:
You can select one of possibilities listed above using options fwhm-constrain and shape-constrain.
Phase is defined by its type and lattice parameters. Command
adds a new phase to set of phases in our model. p is a letter that stands for type of lattice:
%c - cubic phase,
%t - tetragonal phase,
%o - orthorhombic phase,
%h - hexagonal phase.
Monoclinic and triclinic types are not implemented, because author of program had no need to use them. You can also display this list with command:
Lattice parameters, that follow type of phase, can be specified in the same way as lambda and intens in c.wavelength command described above.
The last information, which the program needs, is what planes are visible on diffractogram. It has to be specified using Miller indices (hkl). Command
adds appropriate functions, that represent that plane, to the sum of functions. If %n is omitted, first phase (%0) is assumed. Eg. command c.add %1 (100) will add plane (100) of second defined phase %1. Informations about phases and planes can be viewed with command:
Every plane or the whole phase can be removed from the model using:
The last thing that requires explanation is how initial height, width and shape parameters of peaks are being found. It is based on algorithm used in m.findpeak command. It is also influenced by the same options as m.findpeak. In case of two or more wavelengths used, this algorithm is a little enhanced. To see, how peaks representing plane will look like, use command:
If width is given, it is used instead of value of search-width option.
You can find following script and data used by it in samples directory.
## this script is an example of crystallography/diffraction module usage d.load 'SiC_Zn.dat' # Load data file from current directory (SiC+Zn data) #background was obtained by clicking on plot in "background mode" d.background 20.4823 43.3474 , 28.004 24.3128 , 31.4412 23.6984 , d.background 33.7911 36.0138 , 38.6477 30.5403 , 98.9835 23.6079 , d.b 98.9835 23.6079 , 79.5965 22.1979 , 65.1888 17.968 , 49.9335 17.263 d.range [ 31 : ] # only data with 2theta greater then 31.0 are active #Used wavelength in known and we do not want to fit it. # If you would like to fit wavelength, replace '_' with '~' c.w _1.54051 # Cu Ka1 c.set peak-type = Mod-TCHpV #define phases and initial lattice parameters c.add %c ~4.35 #SiC, %0 c.add %h ~2.66 ~4.92 #Zn, %1 #defining some peaks (order does not matter) c.add %0 (220) # (220) plane of SiC phase c.add %1 (002) # (002) plane of Zn phase c.add %0 (111) c.add %1 (100), %1 (101) # do you remember how comma(,) works? c.a %1 (102) # every command can be shortened with x.x (letter-dot-letter) #after adding a few peaks, it can be a good idea to fit it, to obtain #better approximation of lattice parameters #f.run 10 # the rest of peaks c.add %1 (103), %1(110), %1(112), %1(200), %1(201), %1(202) c.add %0 (311) s.add <x ~0 #zero-shift f.run # running default fitting method (Lev-Mar) #'lambda' that you can see in output window is a parameter used by # fitting method, not wavelength. c.info # see refined lattice parameters # if you want to zoom, the easiest way is to use left button on auxiliary # plot (this plot with difference). To zoom out, use middle button on it. # To see individual peaks, select Show->Peaks from pop-up menu.
[1] Syntax of this command will be changed, other error estimates or measures of goodness-of-fit will be added.
autoplot (o.set) |
background type (d.set) |
cancel-peak-out-of-search (m.set) |
cut-tails (s.set) |
cut-tails-level (s.set) |
default-max-iterations (f.set) |
default-relative-domain-width (s.set) |
distrib-type (f.set) |
estimate-consider-sum (m.set) |
exit-on-warning (o.set) |
fwhm-constrain (c.set) |
lambda-down-factor (f.set) |
lambda-starting-value (f.set) |
lambda-up-factor (f.set) |
max-wssr-evaluations (f.set) |
min-fract-range (f.set) |
minimal BG points distance (d.set) |
move-all (f.set) |
move-multiplier (f.set) |
numerical-d (s.set) |
numerical-d-both-sides (s.set) |
numerical-d-rel-change (s.set) |
output-one-of (f.set) |
peak-type (c.set) |
pseudo-random-seed (f.set) |
recursive-remove (s.set) |
search-width (m.set) |
shake-before (f.set) |
shake-type (f.set) |
shape-constrain (c.set) |
stop-rel-change (f.set) |
verbosity (o.set) |
To simplify equations, following substitution is introduced:
where a1 is a position of peak's center and a2 is a HWHM.
a3 parameters in Pearson VII and Pseudo-Voigt are not related.
Voigt function is a convolution of Gaussian and Lorentzian functions. a0 = heigth, a1 = center, a2 is proportional to the Gaussian width, and a3 is proportional to the ratio of Lorentzian and Gaussian widths. And t is just a variable of integration. Voigt is computed according to R.J.Wells, “ Rapid approximation to the Voigt/Faddeeva function and its derivatives ”, Journal of Quantitative Spectroscopy & Radiative Transfer 62 (1999) 29-48. (See also: http://personalpages.umist.ac.uk/staff/Bob.Wells/voigt.html). Is the approximation exact enough for all possible uses of fityk program?
weighted sum of squared residuals
Command Line Interface
Graphical User Interface
Full Width at Half Maximum
Half Width at Half Maximum
Fityk is free software; you can redistribute and modify it under terms of GNU General Public License, version 2. There is no warranty. GPL is one of most popular licenses, and it is worth reading, if you have not done it before. The program is copyrighted by author, and the license itself by FSF. Text of license is distributed with program in file COPYING.
This manual is written by author of program in DocBook (XML) and converted to other formats. All changes, improvements, fixes of mistakes, etc. are welcome. The fitykhelp.xml file is distributed with program sources, and can be modified with any text editor.
[1] Numerical Recipes in C. http://www.nr.com.
[2] Data Fitting in the Chemical Sciences by the Method of Least Squares . John Wiley & Sons. 1992.
[3] Data Analysis. Springer Verlag. 1999.
[4] PeakFit 4.0 for Windows User's Manual. AISN Software. 1997.
[5] Algorytmy genetyczne + struktury danych = programy ewolucyjne. WNT. 1996.
[6] R. A. Young. The Rietveld Method. Oxford University Press. 1993.
[7] User's Guide to Program DBWS-9807a. 2000.