Multi-FRILLS is a subroutine library designed to be used in programs that perform interactive non-linear least-squares fitting of a parametrized function to one or more sets of data points. It provides a flexible and easy-to-use interface to a least-squares fitting routine with the minimum of programming on the user's part. The simplest use of Multi-FRILLS is to fit a function to a single set of data points, e.g. counts as a function of time in a single detector. More sophisticated programs can fit several data sets simultaneously, with some parameters global to all data sets and others specific to a particular data set.
The user is required to provide:
When subroutine MULTI_FRILLS has been called, the program operates in command mode. Following the prompt '#', the user types commands to
On completion, the fitted parameters are listed along with their estimated errors and a value for chi^2. The program returns to command mode so that the user can plot the fit and choose new initial conditions, e.g. by clearing fixed parameters, before trying again. The results of fits are recorded in an ASCII file and lists of parameters may be stored in files for subsequent printing or input as initial parameters to another Multi-FRILLS session. The plots may also be stored on disk in Postscript or LN03 files for output to a printer, or the calculated function may be written to ASCII files in tab-delimited form suitable for input to PC or Macintosh programs such as Excel, SigmaPlot, proFit or Kaleidagraph.
To optimize the function parameters, Multi-FRILLS uses either its own non-linear least-squares fitting routine which employs a standard quasi-Newton (inverse Hessian) algorithm, or a standard Levenberg-Marquadt algorithm. For a discussion of the general principles, see section 14.4 of "Numerical Recipes: The Art of Scientific Computing" (Cambridge University Press 1986). Using the calculated values passed back from the user-supplied subroutine, a chi^2 function is defined in the following way:
where v and s are the data values and errors, c are the calculated values, and nv and np are the number of data values and fitting parameters respectively. The best fit parameters are those which minimize chi^2. With this definition, the minimum of chi^2 should be close to
Much larger values of chi^2 mean either that the model function is not an accurate representation of the data or that the data errors are underestimated, whereas smaller values mean that the data errors are too large. The fitting method will find the local chi^2 minimum; there may be a global minimum elsewhere in parameter space. However, since the user can plot the calculated function and compare it to the data for each data set before attempting a fit, it should be possible to define reasonable starting parameters for most models. The criteria for chi^2 convergence may be adjusted along with other fitting constants such as the maximum number of iterations. See the description of the ALTER FIT command in the next section for more details. Estimates of the parameter errors are provided by the diagonal elements of the covariance matrix, which assumes that they are normally distributed i.e. uncorrelated. The covariance matrix can be inspected for evidence of such correlations.
Multi-FRILLS v1.3 is based on FRILLS which is described in the report
"FRILLS: An Interactive Least Squares Fitting Package"
R. Osborn, Rutherford Appleton Laboratory Report RAL 91-011 (1991)
At present, Multi-FRILLS can only be installed on VMS machines on which Dec GKS has been installed. We are in the advanced stages of a project to convert the graphics to PGPLOT, which is a freely available library of Fortran graphics routines. It should then be possible to install Multi-FRILLS on platforms other than VMS.
Section 2 describes all the Multi-FRILLS commands in detail; if you are just using a Multi-FRILLS application and not writing your own, this is the only section you need to read. For programmers, a description of how to write Multi-FRILLS applications is given in Section 3 together with a simple example. A Multi-FRILLS application may launched from the data analysis package Open GENIE (or the older GENIEv2), with the data being transferred from one or more GENIE workspaces. A description of how to do this is given in Section 4.
An important feature of this fitting package is the flexibility that the user has in the way a fit is performed. Commands may be given in any order, with the minimum of constraints. For instance, parameters may be fixed before one fit and cleared before the next. At any stage, the function, defined by the current set of parameters, may be plotted on the terminal screen for any of the several data sets. This section describes all the currently available Multi-FRILLS commands for setting up the parameters, plotting the data and functions, and performing the least-squares fits.
The following conventions are followed by Multi-FRILLS:
In this manual, optional parameters are surrounded by square brackets. Values to be chosen by the user are surrounded by angle brackets.
Multi-FRILLS graphics work both on X-Windows screens and on terminals/emulators with both text and Tektronix 4014 graphics screens, e.g. VersaTerm Pro. On the latter, the terminal is automatically switched to the graphics screen after a plot command. However, the user must type [RETURN] after the plot is completed in order to return to text mode.
At any stage in a Multi-FRILLS session, help can be obtained using the HELP command:
Format: H [<command>]
Format: U <n>
Format: S [<n1>,[<n2-n3>,[n4...]]]
Format: F <n1>,[<n2-n3>,[n4...]]
Format: B <n1>,[<n2-n3>,[n4...]]
Format: L <n1>,[<n2-n3>,[n4...]]
Format: I <disk file name>
Format: O P [<disk file name>]
Format: O F [<disk file name>]
Format: R <xmin> <xmax>
Format: M <xmin> <xmax>
Format: C P [<n1>,[<n2-n3>,[n4...]]]
Format: C D [<xmin> <xmax>]
Format: D D [<xmin> <xmax>]
Format: D P [<n1>,[<n2-n3>,[n4...]]]
Format: D C [<xmin> <xmax>]
Format: D F
Format: D M
Format: D L
Format: P [O] D [<xmin> <xmax> [<ymin> <ymax>]]
Format: P [O] C [<xmin> <xmax> [<ymin> <ymax>]]
Format: P [O] F [<xmin> <xmax> [<ymin> <ymax>]]
Format: P [O] R [<xmin> <xmax> [<ymin> <ymax>]]
Format: P [O] [<n1>,[<n2-n3>,[n4...]]]
Format: O C [<disk file name>]
Format: K [<disk file>]
Format: A F
The last parameter illustrates where the quasi-Newton fitting algorithm used as the default in Multi-FRILLS parts company with other common non-linear least-squares methods. When the inverse Hessian method produces an increase in chi^2, Multi-FRILLS halves the most recently calculated parameter steps until a smaller value of chi^2 results. The minimization then continues using the new values. This simple procedure has proved to be very successful in the sorts of problems most frequently solved by Multi-FRILLS. However, the more commonly-used Levenberg-Marquadt algorithm is available as an option.
Format: A O [<ncode>]The possible values of
ncode are:
0 No output to terminal during fit.
1 Output chi^2 to terminal during fit.
2 Output chi^2 and params. to terminal during fit.
10 As for option 0, but, in addition, print fitted points in file.
11 As for option 1, but, in addition, print fitted points in file.
12 As for option 2, but, in addition, print fitted points in file.
Format: A M [<ntype>, [<size> ], <ncolour>]]]
Format: A L [<ntype>, [<thickness> ],]]]
Type [RETURN] to leave them unchanged. The function is calculated at N_PLT (default 1000) x-coordinates equally spaced over the range of the plot. N_PLT should be reduced if the function is computationally intensive to evaluate. If the user sets N_PLT=1, the function is calculated only at those x-coordinates at which there are data and which have not been REMOVEd. If the user sets SMOOTH=1, then from the N_PLT calculated data points a smooth curve is generated at 1000 equally spaced x-coordinates which passes through the calculated data points. This produces a smooth curve even if the function is evaluated at only a few points but may be an inaccurate representation of the true function if the number of nodes is small. By default SMOOTH=0 and the calculated values are joined by straight lines.
Format: A S
N.B. the smoothing options only affect the appearance of the calculated function in plots. chi^2 is only calculated using the data points themselves, so the fitting results are not affected by these options.
Format: A W
Format: A H <ntype>
Format: T <new-title>
Format: J [<DCL command>]
Format: G
Format: E
In this section all the essential initialization required by Multi-FRILLS is described, together with the input and output variables of the fitting function subroutine. An example application and an example session are also provided. The source code for the example program and the GENIE examples in Section 4, is provided with the Multi-FRILLS source code and installation package. Following installation of the package, two logical names should be defined:
If these have not been defined at the system level, please consult your system administrator, or whoever installed the Multi-FRILLS library to find out where the relevant subdirectories are located. You should then add the definitions to your own LOGIN.COM e.g.
$ DEFINE FRILLS_SOURCES DUA0:[SOME.SUBDIRECTORY.SOURCES] $ DEFINE FRILLS_LIBRARY DUA0:[SOME.SUBDIRECTORY.LIBRARY]
The user needs to provide the main program which should initialize the data points, axis labels and the parameter names before calling MULTI_FRILLS. It is also possible to preset initial parameter values and whether they are to be fitted or not. This is useful for those functions in which a number of parameters are never intended to be fitted, e.g. the sample temperature. Note however that these parameters can be cleared at a later stage by the CLEAR PARAMETER command, so the user should be careful in using this feature. The user-supplied subroutine that defines the fitting function, must be explicitly declared as EXTERNAL in the program that calls MULTI_FRILLS e.g.
external GAUSS
The calling program should also contain the line
include 'FRILLS.INC'
which gives access to all the common blocks used by Multi-FRILLS. The user needs to initialize the following variables:
| NPTOT | I*4 | the total number of parameters |
| NAM(ip) | C*20 | the names of the parameters (ip=1:NPTOT) |
| INS(ip) | I*4 | the data sets with which each parameter is associated |
| NS | I*4 | the number of data sets |
For each of the data sets (is=1, NS), the following values must be set:
| ND(is) | I*4 | the number of points in each data set (is=1:NS) |
| X(i,is) | R*8 | the x-coordinates (i=1:ND(is), is=1:NS) |
| Y(i,is) | R*8 | the y-coordinates |
| YSIG(i,is) | R*8 | the error bars |
| XLAB(is), YLAB(is) | C*5 | names used for the x and y axes when DISPLAYing data |
| XCAP(is), YCAP(is) | C*40 | labels used for the x and y axes when PLOTting data |
| TITLE(is) | C*80 | title |
| SUBTITLE(3,is) | C*60 | optional subtitles to be placed under the plots |
Optionally, the user can define each parameter as either global to all data sets or specific to one data set, by filling the array INS(ip) = 0 (global) or INS(ip) = is (specific to data set is). The default is for parameters to be global to all data sets. The array INS is used for display purposes only and is not used in the minimization algorithm.
The Multi-FRILLS command interface is accessed by a call to subroutine MULTI_FRILLS with the fitting function subroutine name as argument e.g.
call multi_frills (GAUSS)
The include file FRILLS_SOURCES:FRILLS.INC gives access to many more variables which allow the ambitious programmer more control over the initial state of Multi-FRILLS. Details are given in the next section. Whilst the use of include files simplifies the writing of Multi-FRILLS applications, it opens up the danger that the user will inadvertently use variable names that are also contained in the Multi-FRILLS common blocks. The simplest way to reduce this possibility is for the user to avoid the IMPLICIT statement and explicitly declare all variables in the calling program. The DEC FORTRAN compiler will print a warning message if any variable has been declared twice and the compiler option FORTRAN/WARNINGS=DECLARATIONS will print a warning if there are any undeclared variables. This option is used in the command files provided with the library. The common blocks are included in the appendix to this manual for reference.
Diagnostic output from the fit is written to unit 41. If this is not explicitly opened by the calling program, it will be sent to SYS$SCRATCH:FRILLS.TXT.
The subroutine defining the fitting function, in following example GAUSS, is expected to calculate the fitting function for only one data set at each call. It has no arguments passed to it directly, but it must contain the statement
include 'FRILLS_SOURCES:FUNCTION.INC'
which gives access to the following input variables:
| P(ip) | R*8 | the current values of the parameters. |
| IS | I*4 | the index of the data set for which calculated values are required. |
| NV | I*4 | the number of points in the data set. |
| X(iv) | R*8 | the x-coordinates of the points (iv=1:NV). |
| X(iv) | R*8 | the x-coordinates of the points (iv=1:NV). |
and the variables to contain the output from the fitting function:
| YCAL(iv) | R*8 | the calculated function at each of the x values (iv=1:NV) |
| NC | I*4 | the number of separate function components that make up YCAL |
| YC(iv, ic) | R*8 | the calculated function components (iv=1:NV; ic=1:NC) |
Filling NC and YC is optional. They permit the individual components of a multi-component function to be plotted separately with the PLOT command. Although the data arrays in FUNCTION.INC only contain one data set, the parameter array contains all the parameters. The subroutine defining the function must distinguish which parameters are required by the data set being processed. We provide an example of how this may be done.
Multi-FRILLS applications can be compiled and linked either on VAX/VMS or Alpha/VMS computers on which DEC GKS has been installed. To compile and link a program to Multi-FRILLS, use the following command:
$ @FRILLS_LIBRARY:FRILLS_LINK <program name(s)> [<library name(s)>]This compiles and links the program both to the Multi-FRILLS object library FRILLS_LIBRARY:MULTI_FRILLS.OLB, and to the DEC Fortran Binding GKS library SYS$LIBRARY:GKSFORBND.OLB. Other object files or libraries can be attached to <program name> by adding the second parameter to the command line. These must explicitly include the /LIB or /OPT qualifiers as they would appear in a normal link statement. Alternatively, the user can edit FRILLS_LIBRARY:FRILLS_LINK.COM to include these libraries automatically. To allow for reasonable numbers of parameters, data sets and the lengths of data sets Multi-FRILLS requires a large amount of computer memory. The program may fail to compile as it might exceed the user's page file quota. Please see your system administrator if this happens.
The Multi-FRILLS library currently has the following defaults:
The number of parameters that are allowed to vary in a fit should be much less than 150 for acceptable running time.
The limits are contained in the file FRILLS_SOURCES:FRILLS.PAR which may be edited. However, the whole Multi-FRILLS library would need to be recompiled before linking the new program. If Multi-FRILLS is installed centrally at your site, consult the system administrator about changing the options. Alternatively, you could download the Multi-FRILLS library in order to have your own customized copy.
The following program is an example of a Multi-FRILLS application to fit a single Gaussian function to a number of data sets. Although the parameters can be different for each data set, Multi-FRILLS allows the user to bind parameters together. For example, if the user believes that the centres of the peak should be the same in each data set, they can all be bound to a common value. In the calling program, two data sets are read in from FRILLS_SOURCES:GAUSS_FIT.DAT, though the rest of the program would accommodate more data sets. The source code is contained in FRILLS_SOURCES:GAUSS_FIT.F. The data sets provided with the library have been produced by defining a Gaussian on top of a sloping backgroun. Statistical fluctuations have been simulated by assuming a Gaussian distribution of errors with variances equal to the data values. For comparison with the fitted parameters, the Gaussian and background parameters were
| Parameter | Dataset 1 | Dataset 2 |
|---|---|---|
| Height | 200 | 100 |
| FWHM | 1.5 | 3 |
| Centre | 5 | 5 |
| Background Constant | 10 | 10 |
| Background Slope | 1 | 1 |
FRILLS_SOURCES:GAUSS_FIT.F
!------------------------------------------------------------------------------------
! Multi-FRILLS - An Interactive Least Squares Fitting Package for Multiple Datasets
!
! Copyright (C) 1996 R. Osborn, T. G. Perring (See MULTI_FRILLS.FOR for details)
!
! Version Date Name Comments
! ======= ==== ==== ========
! 1.0 96/08/23 RO/TGP First release of multiple dataset version of FRILLS
!
!------------------------------------------------------------------------------------
! program GAUSS_FIT
!
! Program to fit single Gaussian lineshape to more than one dataset
!
! Supplied as an example of a simple Multi-FRILLS application
!
program GAUSS_FIT
include 'FRILLS.INC'
external GAUSS
integer*4 i, j, k
!Read in the data
open (unit = 10, file = 'FRILLS_SOURCES:GAUSS_FIT.DAT', status = 'OLD')
ns = 2 !No. of datasets
nd(1) = 41 !No. of points in dataset 1
nd(2) = 41 !No. of points in dataset 2
do i = 1,41
read (10, *) x(i,1), y(i,1), ysig(i,1), y(i,2), ysig(i,2)
x(i,2) = x(i,1)
end do
close (10)
!Define axis labels and titles
xcap(1) = 'Energy Transfer [meV]' !Extended X-caption for graphical output
xcap(2) = xcap(1)
xlab(1) = 'En' !Abbreviated X-label for data columns
xlab(2) = xlab(1)
ycap(1) = 'Neutron Counts' !Extended Y-caption for graphical output
ycap(2) = ycap(1)
ylab(1) = 'Cts' !Abbreviated Y-label for data columns
ylab(2) = ylab(1)
title(1) = 'Sample Dataset 1' !Title for dataset 1
title(2) = 'Sample Dataset 2' !Title for dataset 2
!Loop over no. of datasets defining parameter names and keep count of parameter nos.
i = 0
do j = 1,ns
write (nam(i+1), '(a,i2,a)') 'Set ', j, ' Height'
write (nam(i+2), '(a,i2,a)') 'Set ', j, ' FWHM'
write (nam(i+3), '(a,i2,a)') 'Set ', j, ' Centre'
write (nam(i+4), '(a,i2,a)') 'Set ', j, ' Bkgd Constant'
write (nam(i+5), '(a,i2,a)') 'Set ', j, ' Bkgd Slope'
do k = 1,5
ins(i+k) = j !Parameters are associated with dataset k
inp(i+k) = i + k !All parameters are normally varied
end do
i = i + 5
end do
nptot = i !Total number of parameters
np = i !Number of parameters to be varied
!Open output file to contain diagnostic details of fit
open (unit = 41, file = 'SYS$SCRATCH:GAUSS_FIT.OUT', status = 'NEW')
!Pass control to Multi-FRILLS program
call multi_frills (GAUSS)
!Close output file
close (unit=41)
stop
end
!------------------------------------------------------------------------------------
! subroutine GAUSS
!
! Subroutine defining the Gaussian function to be fitted to the data
! Only one dataset is calculated in each subroutine call identified by 'is'
!
! Function parameters and dependent variables are passed through FUNCTION.INC
subroutine GAUSS
include 'FUNCTION.INC'
integer*4 i, j
real*8 h, w, c, bc, bs, s, const
parameter (const = 2D0*sqrt(2D0*log(2D0)))
!Definition of parameters
i = 5 * (is - 1)
h = p(i+1)
w = p(i+2)
c = p(i+3)
bc = p(i+4)
bs = p(i+5)
!Two components to the fit are defined - the peak and the background respectively
nc = 2
!Loop over data points
s = w / const !Convert FWHM to standard deviation
do i = 1,nv
yc(i,1) = h * exp (-(x(i)-c)**2/(2D0*s**2))
yc(i,2) = bc + bs * x(i)
ycal(i) = yc(i,1) + yc(i,2)
end do
return
end
We reproduce here an annotated "typical" Multi-FRILLS session using the GAUSS-FIT application shown in the previous section (cf the fitted parameters with those used in producing the simulated data). During the Multi-FRILLS session, the user is prompted for commands with '#'.
$ @frills_library:frills_link gauss_fit
$ run gauss_fit Multi-FRILLS version 1.3, Copyright (C) 1996, 1997 R. Osborn, T. G. Perring
The program reads in the data and then enters the Multi-FRILLS interface, indicated by the command prompt #.
# plot dataAll Multi-FRILLS command names and keywords are uniquely identified by their first letter. In this instance it is sufficient to type just P D. The plotting range can be controlled by typing additional parameters on the same line. Try:
# p d 0 10 0 250
# u 2 # p d
# u 1 (Reset to first data set)
# s
Dataset 1
==========
Give the parameter value (with min and max if required)
1 Set 1 Height : 200
2 Set 1 FWHM : 4
3 Set 1 Centre : 4
4 Set 1 Bkgd Constant:
5 Set 1 Bkgd Slope :
# d p
Dataset 1
==========
Parameter set P Pmin Pmax
1 1 Set 1 Height 200.00
2 1 Set 1 FWHM 4.0000
3 1 Set 1 Centre 4.0000
4 1 Set 1 Bkgd Constant 0.00000E+00
5 1 Set 1 Bkgd Slope 0.00000E+00
Note that only the first five parameters are displayed. This is because they are the only ones associated with the first data set. In the output to the DISPLAY PARAMETER command, the first column is the parameter number and the second column is the data set with which it is associated. If any global parameters were defined in this application, these would have 0 in the second column and would be displayed with all the data sets.
# p c
# s 2,3
Give the parameter value (with min and max if required)
2 Set 1 FWHM : 2
3 Set 1 Centre : 4.5
# u 2
# s
Dataset 2
==========
Give the parameter value (with min and max if required)
6 Set 2 Height : 100
7 Set 2 FWHM : 3
8 Set 2 Centre : 4
9 Set 2 Bkgd Constant:
10 Set 2 Bkgd Slope :
# b 8
Give the binding parameter no. and the ratio
Set 2 Centre : 3 1
# d p
Dataset 2
==========
Parameter set P Pmin Pmax
6 2 Set 2 Height 100.00
7 2 Set 2 FWHM 3.0000
8 2 Set 2 Centre 4.5000 Bound to parameter 3
9 2 Set 2 Bkgd Constant 0.00000E+00
10 2 Set 2 Bkgd Slope 0.00000E+00
# f 5,10
# g
Iteration 1
Old Chi**2 = 0.23542033D+02 New Chi**2 = 0.29901154D+01
Iteration 2
Old Chi**2 = 0.29901154D+01 New Chi**2 = 0.12794530D+01
Iteration 3
Old Chi**2 = 0.12794530D+01 New Chi**2 = 0.11186471D+01
Iteration 4
Old Chi**2 = 0.11186471D+01 New Chi**2 = 0.11171174D+01
Iteration 5
Old Chi**2 = 0.11171174D+01 New Chi**2 = 0.11171038D+01
Job completed: Agreement factor converged
Chi ** 2 1.1171
Parameter set P Sigma
1 1 Set 1 Height 190.94 7.7412
2 1 Set 1 FWHM 1.5542 0.47325E-01
3 1 Set 1 Centre 5.0163 0.20700E-01
4 1 Set 1 Bkgd Constant 13.875 0.78519
5 1 Set 1 Bkgd Slope 0.00000E+00 Fixed
6 2 Set 2 Height 103.57 4.2283
7 2 Set 2 FWHM 2.9573 0.11575
8 2 Set 2 Centre 5.0163 Bound to parameter 3
9 2 Set 2 Bkgd Constant 15.225 1.1257
10 2 Set 2 Bkgd Slope 0.00000E+00 Fixed
Do you wish to restore old parameters [N]?
# p f
# c p 5,10
# g
Iteration 1
Old Chi**2 = 0.11477094D+01 New Chi**2 = 0.77440763D+00
Iteration 2
Old Chi**2 = 0.77440763D+00 New Chi**2 = 0.77435676D+00
Job completed: Agreement factor converged
Chi ** 2 0.77436
Parameter set P Sigma
1 1 Set 1 Height 191.10 6.4791
2 1 Set 1 FWHM 1.5393 0.39569E-01
3 1 Set 1 Centre 4.9978 0.17568E-01
4 1 Set 1 Bkgd Constant 10.640 0.99233
5 1 Set 1 Bkgd Slope 0.79141 0.18196
6 2 Set 2 Height 102.95 3.5557
7 2 Set 2 FWHM 2.9139 0.96785E-01
8 2 Set 2 Centre 4.9978 Bound to parameter 3
9 2 Set 2 Bkgd Constant 12.206 1.1879
10 2 Set 2 Bkgd Slope 0.81031 0.19721
Do you wish to restore old parameters [N]?
# u 1 # p f # a l 2 # p o 1 # a l 3 # p o 2 # k gauss1.ps Plot stored in UD1:[OSBORN]GAUSS1.PS;1
# u 2 # a l 1 # p f # a l 2 # p o 1 # a l 3 # p o 2 # k gauss2.ps Plot stored in UD1:[OSBORN]GAUSS2.PS;2 # o f gauss.fit Fitted parameters stored in UD1:[OSBORN]GAUSS.FIT;1 # e
Open GENIE is an interactive data analysis package developed at the ISIS Pulsed Neutron Facility, Rutherford Appleton Laboratory . It is in widespread use at ISIS, and other neutron facilities, as a means of rapid visualization and simple analysis of neutron scattering data. The neutron data are grouped into workspaces which contain the x and y-values, errors, units, titles etc of an individual spectrum. We have written a number of subroutines to ease the transfer of these workspaces into the Multi-FRILLS applications, automatically copying the data, graph labels and titles into the Multi-FRILLS common blocks.
Open GENIE supercedes GENIE v2, which would only operate on VMS machines and had a more limited command language. However, the old version of GENIE is still supported by Multi-FRILLS and will also be described here. Source code for both versions are stored in the examples subdirectory i.e. if FRILLS_LIBRARY points to DUA0:[MULTI_FRILLS], the examples will be in either DUA0:[MULTI_FRILLS.EXAMPLES.OPENGENIE] or DUA0:[MULTI_FRILLS.EXAMPLES.GENIE].
In the following sections, we discuss how to run, write and link GENIE applications. Finally, we give an example of a program that fits several GENIE workspaces at once. The source code is in the example file MULTI_RESFIT.F. The program models the peak shape in all workspaces as the convolution of a Gaussian with an exponential decay with a sloping background. A number of the functions used to define the final peak shape are included with other common peak functions in the file FRILLS_SOURCES:PEAK_FUNCTIONS.F which is also included with the Multi-FRILLS package. Other examples of GENIE applications also included are:
To run a Multi-FRILLS application from within Open GENIE, the user should first load the GCL command file:
>> load "frills_library:frills_genie.gcl"
In order to run the program, the user types
>> <result> = frills_genie("<program name>", <workspace>)
where
The operation of the Multi-FRILLS application is identical to a standalone application once it has been started.
To run a GENIE application from within GENIE v2, the user should run a GENIE command file:
>> @frills_library:frills_geniewhich will prompt the user for
Generally, Multi-FRILLS programs for GENIE can act on any number of workspaces up to the maximum number for which the program was designed. The operation of the Multi-FRILLS application is identical to a standalone application once it has been started.
It is possible to write GENIE applications that explicitly act on one workspace only. These should be invoked directly by the GENIE TRANSFORM command:
>> TRANSFORM W<m> <program-name> W<n>
Open GENIE applications are written as subroutines which are then linked as shareable objects and loaded into an Open GENIE session. Several general purpose subroutines have been provided so that an Open GENIE application is as simple to write as a stand-alone Multi-FRILLS application. In fact, it is normally simpler, because the data, graph captions and titles are read in automatically from the GENIE workspaces.
In the following, we shall use MULTI_RESFIT.F as an example. The first lines of the application should contain the following :
subroutine MULTI_RESFIT (pars_get, pars_put)
include 'FRILLS.INC'
include 'GENIE_MODULES.INC'
external RESFUN, pars_get, pars_put
PARS_GET and PARS_PUT are used by Open GENIE to pass the workspace information to the Multi-FRILLS application. RESFUN is the name of the subroutine containing the fitting function.
It is advisable to check that the application is consistent with the version of Open GENIE being run.
if (module_version_ok (genie_major, genie_minor) .eq. 0) then
call module_error ("MULTI_RESFIT", "Error: Open GENIE Version mismatch",
+ "Please recompile this module")
return
end if
The simplest way of inputting GENIE workspaces into the application is to replace the call to MULTI_FRILLS with one to MULTI_FRILLS_OPENGENIE
call multi_frills_opengenie (RESFUN, pars_get, pars_put)
This transfers data (assumed to be in histogram mode), including titles and graph captions, from PARS_GET to the Multi-FRILLS common blocks, and then calls MULTI_FRILLS itself. After completion of the Multi-FRILLS session it transfers the function, calculated with the latest set of parameters, into GENIE workspaces. The user need only provide the names and number of parameters, i.e.
| NPTOT | the total number of parameters |
| NAM(ip) | the names of the parameters (ip=1:NPTOT) |
and optionally fill the array elements INS(ip) = 0 (global) or INS(ip) = is (specific to data set is). The subroutine that calculates the fitting function should be written exactly as for standalone Multi-FRILLS applications. For those users who wish to manipulate or inspect the data before MULTI_FRILLS is called, two routines, MULTI_FRILLS_OPENGENIE_IN and MULTI_FRILLS_OPENGENIE_OUT, are provided. In this case the calling sequence would be
call multi_frills_opengenie_in (pars_get)
.
.
call multi_frills (FUN)
.
.
call multi_frills_opengenie_out (FUN, pars_put)
where FUN is the name of the fitting function subroutine. One reason for using this sequence is if the number of workspaces being fitted is not known in advance. The call to MULTI_FRILLS_OPENGENIE_IN determines the value of ns after which the parameters can be initialized. This is illustrated in the example program. Otherwise the GENIE program is written exactly as with MULTI_FRILLS_OPENGENIE.
Since Open GENIE modules are subroutines rather than standalone programs, they should end with
return
end
rather than a STOP statement.
IMPORTANT: Don't ever put a STOP statement in an Open GENIE application. It will terminate the whole Open GENIE session, not just the Multi-FRILLS application. One unfortunate side effect of this method of invoking external procedures is that a fatal error in the procedure could prove fatal to the whole Open GENIE session. Extra care must therefore be taken in writing the fitting functions so that they do not crash. If in doubt, save essential data into an Open GENIE file before running FRILLS_GENIE.
GENIE v2 applications are very similar to regular Multi-FRILLS programs. Firstly, they are written as programs, rather than subroutines, which are spawned from a GENIE v2 session. In other words, the first lines would be :
program MULTI_RESFIT
include 'FRILLS.INC'
external RESFUN
where RESFUN is the subroutine containing the fitting function. The simplest way of inputting GENIE workspaces into the application is to replace the call to MULTI_FRILLS with one to MULTI_FRILLS_GENIE
call multi_frills_genie (RESFUN)
This reads data (assumed to be in histogram mode), including titles and graph captions, from a file written when the command file FRILLS_GENIE is invoked within GENIE. The data is read into the MULTI_FRILLS common blocks before MULTI_FRILLS itself is called within MULTI_FRILLS_GENIE. After completion of the Multi-FRILLS session it transfers the function, calculated with the latest set of parameters, into GENIE workspaces by FRILLS_GENIE. The user need only provide the names and number of parameters, i.e.
| NPTOT | the total number of parameters |
| NAM(ip) | the names of the parameters (ip=1:NPTOT) |
and optionally fill the array elements INS(ip) = 0 (global) or INS(ip) = is (specific to data set is). The subroutine that calculates the fitting function should be written exactly as for standalone Multi-FRILLS applications. For those users who wish to manipulate or inspect the data before MULTI_FRILLS is called, two routines, MULTI_FRILLS_GENIE_IN and MULTI_FRILLS_GENIE_OUT, are provided. In this case the calling sequence would be
call multi_frills_genie_in
.
.
call multi_frills (FUN)
.
.
call multi_frills_genie_out (FUN)
where FUN is the name of the fitting function subroutine. One reason for using this sequence is if the number of workspaces being fitted is not known in advance. The call to MULTI_FRILLS_GENIE_IN determines the value of ns after which the parameters can be initialized. Otherwise the GENIE program is written exactly as with MULTI_FRILLS_GENIE.
Multi-FRILLS applications that act on one workspace only, i.e. which are invoked with the GENIE TRANSFORM command rather than the call to FRILLS_GENIE, can be written with SINGLE_FRILLS_GENIE replacing the call to MULTI_FRILLS. Alternatively, if the user wishes to manipulate the data before entering Multi-FRILLS, the following calling sequence should be used:
call single_frills_genie_in
.
.
call multi_frills (fun)
.
.
call single_frills_genie_out (fun)
An Open GENIE application is written as a subroutine that needs to be linked as a shareable object module to Open GENIE. This is achieved by compiling and linking the source code with the command file OPENGENIE_LINK.COM e.g.
$ @FRILLS_LIBRARY:OPENGENIE_LINK <program name> [<library names>].The module can then be loaded with the MODULE/LOAD command. However, this is done automatically by FRILLS_GENIE.GCL, the procedure invoked when running the Multi-FRILLS application.
The linking procedure for GENIE v2 applications is similar to that for standalone Multi-FRILLS programs except that the GENIE library needs to be included by adding GENIE_LIB/LIB after the Multi-FRILLS library. The command file, GENIE_LINK.COM, has been written to compile and link GENIE applications.
$ @FRILLS_LIBRARY:GENIE_LINK <program name> [<library name(s)>]As with FRILLS_LINK.COM, GENIE_LINK.COM may be edited to add more object libraries to the link statement.
The following program fits a form of the resolution function of pulsed neutron chopper spectrometers that works reasonably well in many cases. It has been written to accommodate multiple data sets although there are no global parameters. One use of such a multiple data set version would be to fit the elastic peaks from a number of different vanadium runs with the resolution parameters bound to a common value. The program is invoked by invoking the procedure FRILLS_LIBRARY:FRILLS_GENIE.GCL from within an open GENIE session. The main program calls MULTI_FRILLS_OPENGENIE_IN first in order to determine the number of data sets for which parameters need to be initialized. It then calls MULTI_FRILLS itself to transfer control to the Multi-FRILLS application. Finally, MULTI_FRILLS_OPENGENIE_OUT is called to transfer the results of the fit back to GENIE before exiting.
MULTI_RESFIT.F
!------------------------------------------------------------------------------------
! Multi-FRILLS - An Interactive Least Squares Fitting Package for Multiple Datasets
!
! Copyright (C) 1996 R. Osborn, T. G. Perring (See MULTI_FRILLS.FOR for details)
!
! Version Date Name Comments
! ======= ==== ==== ========
! 1.0 96/08/23 RO/TGP First release of multiple dataset version of FRILLS
! 1.1 97/08/04 RO Removed VMS logical name from include statement
! 1.2 97/09/08 TGP Addition of Levenberg-Marquadt option
! 1.3 99/01/06 RO/TGP First version compatible with Open GENIE
!
!------------------------------------------------------------------------------------
! subroutine MULTI_RESFIT
!
! Fitting of chopper spectrometer resolution function to multiple GENIE workspaces
!
subroutine MULTI_RESFIT (pars_get, pars_put)
include 'FRILLS.INC'
include 'GENIE_MODULES.INC'
external RESFUN, pars_get, pars_put
integer*4 i, j, k
!Check that version of open Genie is OK
if (module_version_ok (genie_major, genie_minor) .eq. 0) then
call module_error ("MULTI_RESFIT", "Error: Open GENIE Version mismatch",
+ "Please recompile this module")
return
end if
! Give a welcome message
write (*,*) '#'
write (*,*) '# MULTI_RESFIT - Chopper Spectrometer Resolution Fitting Program'
write (*,*) '#'
write (*,*) '# Type H for a list of the available commands'
write (*,*) '#'
write (*,*) '# Details of fits are stored in SYS$SCRATCH:MULTI_RESFIT.TXT'
write (*,*) '#'
! Read in GENIE workspaces
call multi_frills_opengenie_in (pars_get)
! Initialize fitting parameters for each dataset
i = 0
do j = 1,ns
write (nam(i+1), '(a,i2,a)') 'Set ', j, ' Bkgd Constant'
write (nam(i+2), '(a,i2,a)') 'Set ', j, ' Bkgd Slope'
write (nam(i+3), '(a,i2,a)') 'Set ', j, ' Intensity'
write (nam(i+4), '(a,i2,a)') 'Set ', j, ' Gamma'
write (nam(i+5), '(a,i2,a)') 'Set ', j, ' Sigma'
write (nam(i+6), '(a,i2,a)') 'Set ', j, ' Centre'
do k = i+1, i+6
ins(k) = j
end do
i = i + 6
end do
nptot = i
! Output is directed to SYS$SCRATCH:MULTI_RESFIT.TXT
open (unit = 41, file = 'SYS$SCRATCH:MULTI_RESFIT.TXT', status = 'NEW')
! Call MULTI_FRILLS_GENIE to initialize data using GENIE workspace and call Multi-FRILLS
call multi_frills (RESFUN)
! Write out GENIE workspaces
call multi_frills_opengenie_out(RESFUN, pars_put)
! Close the file containing the fitting details
close (unit=41)
return
end
!------------------------------------------------------------------------------------
! RESFUN
!
! Subroutine to calculate fitting function
!
! The function contains two components
! 1) Decaying exponential convolved with a Gaussian defining the peak shape
! 2) Sloping background
!
subroutine RESFUN
include 'FUNCTION.INC'
real*8 bc, bs, int, gamma, sigma, centre, EXP_DECAY
integer*4 i
! Definition of parameters
i = (is-1) * 6
bc = p(i+1) !Background constant
bs = p(i+2) !Background slope
int = p(i+3) !Peak intensity
gamma = p(i+4) !Peak decay constant
sigma = p(i+5) !Peak Gaussian sigma
centre = p(i+6) !Peak centre (1st moment)
!Loop over data points defining peak function and two components
nc = 2
do i = 1,nv
yc(i,1) = int * EXP_DECAY (x(i)-centre, gamma, sigma)
yc(i,2) = bc + bs*x(i)
ycal(i) = yc(i,1) + yc(i,2)
end do
return
end
!------------------------------------------------------------------------------------
! EXP_DECAY
!
! Returns value of decaying exponential convolved with a Gaussian
!
real*8 function EXP_DECAY (eps, gamma, sigma)
real*8 eps, gamma, sigma, GAUSS, x, ERFC, pi
parameter (pi = 3.1415926535897932)
if (sigma .le. 0.0 .and. gamma .le. 0.0) then
EXP_DECAY = 0.0
else if (gamma .le. 0.0) then
EXP_DECAY = GAUSS (eps, sigma)
else if (sigma .le. 0.0) then
x = eps - gamma
if (x .ge. 0.0) then
EXP_DECAY = 0.0
else if (x .lt. 50.0D0*gamma) then
EXP_DECAY = exp (x / gamma) / gamma
else
EXP_DECAY = 0.0
end if
else
x = eps - gamma
if (x .lt. 5.0*sigma) then
EXP_DECAY = exp(0.5D0*sigma**2/gamma**2 + x/gamma) * ERFC((sigma/gamma + x/sigma) / sqrt(2D0)) / (2D0*gamma)
else
EXP_DECAY = 0.0
end if
end if
return
end
!------------------------------------------------------------------------------------
! GAUSS
!
! Returns value of normalized Gaussian centred at 0 with standard deviation SIGMA
!
real*8 function GAUSS (x, sigma)
real*8 x, sigma, pi
parameter (pi = 3.1415926535897932)
GAUSS = exp(-x**2/(2D0*sigma**2)) / (sigma * sqrt(2D0*pi))
return
end
!------------------------------------------------------------------------------------
! ERFC
!
! Complementary Error Function
! Algorithm from Abramovitz and Stegun p.299
! Accurate to 3E-07 for both +ve and -ve x.
!
real*8 function ERFC (x)
real*8 a1, a2, a3, a4, a5, a6, x, z
data a1, a2, a3 /0.0705230784, 0.0422820123, 0.0092705272/
data a4, a5, a6 /0.0001520143, 0.0002765672, 0.0000430638/
z = abs(x)
ERFC = 1.0 + a1*z + a2*z**2 + a3*z**3 + a4*z**4 + a5*z**5 + a6*z**6
if (ERFC .ge. 10.0) then
ERFC = 0.0
else
ERFC = 1.0 / ERFC**16
end if
if (x .lt. 0.0) ERFC = 2.0 - ERFC
return
end
Comments to: Ray Osborn <ROsborn@anl.gov>
Revised: Wednesday, October 27, 1999
|
|