Simfit: simulation, statistical analysis, curve fitting and graph plotting.

User defined models

The SIMFIT package has a comprehensive library of models for simulating or fitting, which will cover most situations. However, you may have to construct your own user-defined models in cases involving procedures like the following ones.

Systems of nonlinear differential equations and Jacobians
Sets of linked or independent nonlinear equations
Evaluation of special functions
Probability integrals and inverses
Integrating one or more functions of one or several variables
Locating roots of nonlinear equations in one or several variables
Constrained nonlinear optimization using partial derivatives
Evaluating or fitting convolution integrals
Models with swap over points or discontinuities

Introduction
Global storage and sub-models
Special functions
Operations with vectors
Integer and logical functions
Using standard expressions
Summary of commands
Examples

1. Introduction

Any serious method for supplying models as ASCII text files should use a readily understood language, it must accept all usual maths operations, including trigonometric and hyperbolic functions, and it should allow the use of special functions and calls to other library (DLL) routines. The SIMFIT curve fitting and simulation programs can use such models, as long as the ASCII files are formatted correctly and the syntax has been checked using program usermod.

Using standard expressions (New feature at Version 7.1.6)

Standard mathematical expressions can be used at any time in the model defining section of a user-supplied file as long as they are contained in a begin{expression} ... end{expression} structure. For example the code

     begin{expression}
     f(1) = p(1) + p(2)x + p(3)x^2 + p(4)x^3
     end{expression}

would define model number 1, i.e. f(1), as a cubic using four variable parameters p(1), p(2), p(3), and p(4), and an independent variable x. When a model file is analysed by SIMFIT it accepts the input to create a virtual stack in reverse Polish, that is, last in first out, or post fix notation. However, when a begin{expression} token is encountered it transforms the equations encountered from standard mathematical equations into reverse Polish until an end{expression} token closes the environment. So, for most users, there would never be any need to learn reverse Polish and models can simply be created in usual mathematical notation using program USERMOD.

Note that models involving more advanced operations such as quadrature, root finding, Chebyshev expansions, convolution integrals, vector algebra, or conditional execution, etc. may require combining standard notation with explicit reverse Polish.

Using reverse Polish (i.e., postfix, or calculator) notation

When developing your own model, you must understand that a reverse Polish (last-in-first-out, or post-fix) procedure works by making a stack such that, at the end of the stack operations, the stack has only one element, which is the value of the model. The test files have comments (to the right of the commands, or following the % characters) to help you understand the way models of any degree of complexity can be constructed.

The functions that can be used (see usermod1.tf1, etc.)

In addition to all the standard mathematical functions you can use functions such as erfc, the gamma function, normal integrals, etc. The list has now been extended to include calls to numerical routines, loops, conditionals, iterations, etc.

The SIMFIT method is unique in that it allows the user to optimise the execution stack by using duplicate, pop and so on, and by allowing many functions to be defined at the same time. However, since few people program calculators or write in PostScript nowadays, it may seem hard to realise that anything can be programmed without using parentheses if you use post-fix. Many students with no programming experience have been able to write files for systems of differential equations after a few hours practise so it can't be that hard. Start with usermod1.tf1, which is a simple line y = mx + c. Read the file then input it into program usermod to evaluate, plot, find zeros, estimate areas, etc. and it will all suddenly become clear. Just do it.

Rules for supplying a user defined function for simulation or fitting

Please observe the use of the special symbol % in this file. The symbol % starting a line is an escape sequence to indicate a change in the meaning of the input stream, e.g. from text to parameters, from parameters to model instructions, from model to Jacobian, etc. Characters occurring more than eight places after the first non-blank character are interpreted as comments and text here is ignored when the model is parsed. The % symbol MUST be used to indicate:-

       i)   start of the file
       ii)  start of the model parameters
       iii) start of the model equations
       iv)  end of model equations (start of Jacobian with diff. eqns.)

The file you supply must have EXACTLY the format now described.

The file must start with a % symbol indicating where text starts The next lines must be the name/details you choose for the model. This would normally be at least 4 not greater than 24 lines. This text is only to identify the model and is not used by simfit. The end of this section is marked by a % symbol. The next three lines define the type of model.
The first of these lines must indicate the number of equations in the model, e.g. 1 equation, 2 equations, 3 equations, etc.
The next must indicate the number of independent variables as in:- 1 variable, 2 variables, 3 variables, etc. or else it could be differential equation to indicate that the model is one or a set of ordinary differential equations with one independent variable.
The next line must define the number of parameters in the model.
With differential equations, the last parameters are reserved to set the values for the integration constants y0(i), which can be either estimated or fixed as required. For example, if there are n equations and m parameters are declared in the file, only m-n can be actually used in the model, since y0(i) = p(m-n+i) for i = 1, 2, ..., n.
Lines are broken up into tokens by spaces.
Only the first token in each line matters after the model starts.
Comments begin with % and are added just to explain what's going on.
Usually the comments beginning with a % can be omitted.
Critical lines starting with % must be present as explained above. The model operations then follow, one per line until the next line starting with a % character indicates the end of the model.
Numbers can be in any format, e.g. 2, 1.234, 1.234E-6, 1.234E6
The symbol f(i) indicates that model equation i is evaluted at this point.
Differential equations can define the Jacobian after defining the model. If there are n differential equations of the form
```
    dy(i)/dx = f(i)(x, y(1), y(2), ..., y(n))
    
```
then the symbol y(i) is used to put y(i) on the stack and there must be a n by n matrix defined in the following way. The element J(a,b) is indicated by putting j(n*(b-1) + a) on the stack. That is the columns are filled up first. For instance with 3 equations you would have a Jacobian J(i,j) = df(i)/dy(j) defined by the sequence:
```
    J(1,1) = j(1),  J(1,2) = j(4),  J(3,1) = j(7)
    J(2,1) = j(2),  J(2,2) = j(5),  J(3,2) = j(8)
    J(3,1) = j(3),  J(3,2) = j(6),  J(3,3) = j(9)
```

Permitted operations and the effects produced

 x             : stack -> stack, x
 y             : stack -> stack, y
 z             : stack -> stack, z
 add           : stack, a, b -> stack, (a + b)
 subtract      : stack, a, b -> stack, (a - b)
 multiply      : stack, a, b -> stack, (a*b)
 divide        : stack, a, b -> stack, (a/b)
 p(i)          : stack -> stack, p(i)         ...  i can be 1, 2, 3, etc
 f(i)          : stack, a -> stack  ...evaluate model since now f(i) = a
 power         : stack, a, b -> stack, (a^b)
 squareroot    : stack, a -> stack, sqrt(a)
 exponential   : stack, a -> stack, exp(a)
 tentothepower : stack, a -> stack, 10^a
 ln (or log)   : stack, a -> stack, ln(a)
 log10         : stack, a -> stack, log(a) (to base ten)
 pi            : stack -> stack, 3.1415927
 sine          : stack, a -> stack, sin(a)       ... radians not degrees
 cosine        : stack, a -> stack, cos(a)       ... radians not degrees
 tangent       : stack, a -> stack, tan(a)       ... radians not degrees
 arcsine       : stack, a -> stack, arcsin(a)    ... radians not degrees
 arccosine     : stack, a -> stack, arccos(a)    ... radians not degrees
 arctangent    : stack, a -> stack, arctan(a)    ... radians not degrees
 sinh          : stack, a -> stack, sinh(a)
 cosh          : stack, a -> stack, cosh(a)
 tanh          : stack, a -> stack, tanh(a)
 exchange      : stack, a, b -> stack, b, a
 duplicate     : stack, a -> stack, a, a
 pop           : stack, a, b -> stack, a
 absolutevalue : stack, a -> stack, abs(a)
 negative      : stack, a -> stack , -a
 minimum       : stack, a, b -> stack, min(a,b)
 maximum       : stack, a, b -> stack, max(a,b)
 gammafunction : stack, a -> stack, gamma(a)
 lgamma        : stack, a -> stack, ln(gamma(a))
 normalcdf     : stack, a -> stack, phi(a)  integral from -infinity to a
 erfc          : stack, a -> stack, erfc(a)
 y(i)          : stack    -> stack, y(i)                Only diff. eqns.
 j(i)          : stack, a -> stack  J(i-(i/n), (i/n)+1) Only diff. eqns.
 ***           : stack -> stack, ***          ...  *** can be any number

Error handling

As the stack is evaluated, action is taken to avoid underflow, overflow and forbidden operations, like 1/x as x tends to zero or taking the log or square root of a negative number etc.

2. Global storage and sub-models

This document describes enhancements to the Simfit user-supplied model options at Version 5.4 release 4.037. It is a supplement to w_readme.f5 (which describes the basic techniques) and will be of interest to users who want to define subsidiary models for function evaluation, composition of functions, root finding, adaptive quadrature, constrained optimisation, conditional branching, etc. to be called as independent procedures from within a main model. Note that Simfit creates a separate virtual stack for each model and procedure, and so the only communication between the stacks is the passing of arguments, retrieval of results and sharing of storage and parameter values if requested.

 a) The command put
 b) The command get
 c) The command get3
 d) The command epsabs
 e) The command epsrel
 f) The command blim
 g) The command tlim
 h) The command putpar
 i) The command value
 j) The commands quad and convolute
 k) The command root
 l) The command value3
 m) The command order
 n) The command middle
 o) Syntax for subsidiary models
 p) Rules for subsidiary models
 q) Nesting subsidiary models
 r) IFAIL values for D01AJF, D01EAF and C05AZF
 s) Examples of subsidiary models
 t) Synonyms and abbreviations

a) The command put(i)

This command is used to store the results of a sub-calculation. The top stack element is transferred into storage location i and is available globally to all sub-models called by the main model. The stack length is decreased by one, i.e. the element transferred by put(i), and any existing storage element in location i is overwritten. The main model is checked for use of put(k) with no corresponding get(k) unless the main model calls sub-models. Subsidiary models are not checked, of course.

b) The command get(j)

This command is used to get the results from a sub-calculation. The storage element previously placed in location j is copied onto the top of the stack. The stack length is increased by one, i.e. the element copied by get(j), and storage location j is unchanged. The main model is checked for get(k) with no corresponding put(k). Subsidiary models are not checked, of course.

c) The command get3(i,j,k)

This command is practically equivalent to if ... elseif ... else. If the top stack element is negative get(i) is invoked, if it is zero (to machine precision) get(j) is invoked and, if it is positive get(k) is invoked. The stack length is unchanged, as the top element is replaced by storage location i, j or k. Storage locations i, j and k are unchanged. The main model is checked for corresponding puts. Subsidiary models are not checked, of course. The command order is provided to increase the versatility of the command get3(.,.,.).

d) The command epsabs

This command sets the absolute error for subsequent calculations. The default value of espabs (1.0e-3, i.e. absolute error tolerance) is replaced by the top stack element. The stack length is reduced by one, i.e. the value used to re-define epsabs. If epsabs is used inside a subsidiary model, the effects are confined locally to the submodel.

e) The command epsrel

This command sets the relative error for subsequent calculations. The default value of esprel (1.0e-6, i.e. relative error tolerance) is replaced by the top stack element. The stack length is reduced by one, i.e. the value used to re-define epsrel. If epsrel is used inside a subsidiary model, the effects are confined locally to the submodel.

f) The command blim(i)

This command sets the lower-limit vector for subsequent calculations. The top stack element is used for blim(i), the bottom limit for root finding or lower limit for quadrature. The stack length is reduced by one, i.e. the value used to define blim(i). If blim is used inside a subsidiary model, the effects are confined locally to the submodel.

g) The command tlim(i)

This command sets the upper-limit vector for subsequent calculations. The top stack element is used for tlim(i), the top limit for root finding or upper limit for quadrature. The stack length is reduced by one, i.e. the value used to define tlim(i). If tlim is used inside a subsidiary model, the effects are confined locally to the submodel.

h) The command putpar

This communicates parameters p(i) from the main model to a sub-model. It must be used to transfer the current parameter values into global storage locations if it is wished to use them in a subsidiary model. Unless the command putpar is used in a main model, the sub-models have no access to parameter values to enable the command p(i) to add parameter i to the sub-model stack. The stack length is unchanged. Note that the storage locations for putpar are initialised to zero so, if you do not use putpar at the start of the main model, calls to p(i) in subsequent sub-models will not lead to a crash, they will simply use p(i) = 0. The command putpar cannot be used in a subsidiary model, of course. Note that putpar should be only used when absolutely necessary as it slows down model evaluation.

i) The command value(i)

This command is used to evaluate a subsidiary model. It uses the current values for independent variables to evaluate subsidiary model i. The stack length is increased by one, as the value returned by function evaluation is added to the top of the stack. The command putpar must be used before value(i) if it wished to use the main parameters in subsidiary model number i.
Advice:
It is important to make sure that a subsidiary model is correct, by testing it as a main model, if possible, before using it as a subsidiary model. You must be careful that the independent variables passed to the submodel for evaluation are the ones intended. Of course, value(i) can call sub-models which themselves can call root, and/or quad. At whatever level value(i) is called, the independent variables will be the those at that level of nesting, i.e. top-level from a main model and dummy variables in a sub-model. Expert users will immediately realise that there is nothing to prevent them designing a sub-model that is evaluated without using any parameters or arguments, or perhaps just using gets and/or puts, which gives great scope for ingenuity in defining special functions.

j) The command quad(i)

This command is used to estimate an integral by adaptive quadrature. It uses the epsabs, epsrel, blim and tlim values to integrate model i and place the return value on the stack. The values assigned to the blim and tlim arrays are the limits for the integration. If the model i requires j independent variables, then j blim and tlim values must be set before quad(i) is used. The length of the stack increases by one, the return value placed on the top of the stack. The command putpar must be used before quad(i) if it is wished to use the main parameters in the subsidiary model number i.
Advice:
Adaptive quadrature cannot be expected to work correctly if the range values, e.g. the extreme upper tail of a decaying exponential. The routines used (D01AJF and D01EAF) cannot really handle infinite ranges, but excellent results can be obtained using commonsense extreme limits, e.g. several relaxation times for a decaying exponential, rather than an arbitrarily large number. With difficult problems it will probably be necessary to increase epsrel and epsabs.

The special case of convolution integrals = f*g

These are done using the command convolute(i,j) which has the following special features.

i) convolute(i,j) calculates the integral from A = blim(1) to B = tlim(1) for submodels i (i.e. f(.)) and j (i.e. g(.)) as in
f*g = integral from A to B of f(u)g(B - u)du
ii) Any sub-models can be used for i and j as long as i and j are less than or equal to 10
iii) The lower limit for integration must be blim(1) whatever i or j are
iv) The top limit for integration must be tlim(1) whatever i or j are
v) Extra parameters may need to be set to normalise the integral

The file convolve.mod is an example model illustrating these principles.

k) The command root(i)

This command is used to estimate a zero of a sub-model iteratively. It uses the epsabs, epsrel, blim and tlim values to find a root for model i and places the return value on the stack. The values assigned to blim(1) and tlim(1) are the limits for root location. The length of the stack increases by one, the root value placed on the top of the stack. The command putpar must be used before root(i) if it wished to use the main parameters in the subsidiary model i.

Advice:

The limits A = blim(1) and B = tlim(1) are used as starting estimates to bracket the root. If f(A)*f(B) > 0 then the range (A,B) is expanded by up to ten orders of magnitude (without changing blim(1) or tlim(1)) until f(A)*f(B) < 0. If this or any other failure occurs, the root is returned as zero. Note that A and B will not change sign, so you can search for, say, just positive roots. If this is too restrictive, make sure blim(1)*tlim(1) < 0. C05AZF is used and, with difficult problems, it will probably be necessary to increase epsrel.

l) The command value3(i,j,k)

This is a very powerful command which is capable of many applications of the form: if ... elseif ... else. If the top stack value is negative value(i) is called, if it is zero (to machine precision), value(j) is called, and if it is positive value(k) is called. It relies on the presence of correctly formatted sub-models i, j and k of course, but the values returned by sub-models i, k and k are arbitrary, as almost any code can be employed in models i, j and k. The top value of the stack is replaced by the value returned by the appropriate sub-model. The command order is provided to increase the scope for using the command value3(.,.,.).

m) The command order

This is a mechanism for putting -1, 0 or 1 on the stack to control the effects of the commands get3(.,.,.) and value3(.,.,.). What happens is that, before the command order is used, there must be three numbers on the stack in the order ...,a,w,b where, typically, a and b would be limits, such that a < b, and w would either be below these limits, between these limits, or above these limits. If w =< a then the stack ...,a,w,b will become ...,-1, but if a < w =< b then ...,a,w,b will become ...,0, while, if w > b, then ...,a,w,b will become ...,1. For example, the code

0
x
4
order
value3(1,2,3)
f(1)

is used in the test file updownup.mod to define a model that has two swap-over points where the model definition changes, i.e. sub-model 1 for x =< 0, sub-model 2 for 0 < x =< 4, and sub-model 3 for x > 4. Note that the command order takes three arguments off the stack and replaces them with one argument, so the stack length is decreased by two.

n) The command middle

There must be three numbers on the stack ...,a,w,b as for the command order, with a < b. If w < a the stack will become ...,a, if w > b then the stack will become ...,b, otherwise, for a < w =< b the stack will become ...,w. This command is used to prevent underflow, overflow or calculations that would lead to singularities. For example, the code

0
x
1
middle

would leave a number w on the stack, where 0 =< w =< 1, and w = x only if 0 =< x =< 1. Note that the command order takes three arguments off the stack and replaces them with one argument, so the stack length is decreased by two.

o) Syntax for subsidiary models

The format for defining sub-models is very strict and must be used exactly as now described. Suppose you want to use n independent equations. Then n user files are developed and, when they have been tested, they are pasted in order directly after the end of the main model, each surrounded by a begin{model(i)} and end{model(i)} command. So, if you want to use a particular model as a sub-model, you first of all develop it using program usermod then, when it is satisfactory, just add it to the main model. However, note that sub-models are subject to several restrictions.

p) Rules for subsidiary models

Sub-model files must be placed in numerically increasing order at the end of the main model file. Model parsing is abandoned if a sub-model occurs out of sequence.
There must be no spaces or non-model lines between the main model and the subsidiary models, or between any subsequent sub-models.
Sub-models cannot be differential equations.
Sub-models of the type being described must define just one equation.
Sub-models are not tested for consistent put and get commands, since puts might be defined in the main model, etc.
Sub-models cannot use putpar, since putpar only has meaning in a main model.
Sub-models can use the commands value(i), root(j) and quad(k), but it is up to users to make sure that all calls are consistent.
When the command value(i) is used, the arguments passed to the sub-model for evaluation are the independent variables at the level at which the command is used. For instance if the main model uses value(i) then value(i) will be evaluated at the x, y, z, etc. of the main model, but with model(i) being used for evaluation. Note that y(k) must be used for functions with more than three independent variables, i.e. when x, y and z no longer suffice. It is clear that if a model uses value(i), then the number of independent variables in that model must be equal to or greater than the number of independent variables in sub-model(i).
When the commands root(i) and quad(j) are used, the independent variables in the sub-model numbers i and j are always dummy variables.
When developing models and subsidiary models independently you may get error messages about x not being used, or a get without a corresponding put. Often these can be suppressed by using a pop until the model is developed. For instance x followed by pop will silence any messages about x not being used, etc.

q) Nesting subsidiary models

The subsidiary models can be considered to be independent except when there is a clash that would lead to recursion. For instance, value(1) can call model(1) which uses root(2) to find a root of model(2), which calls quad(3) to integrate model(3). However, at no stage can there be simultaneous use of the same model as value(k), and/or quad(k), and/or root(k). The same subsidiary model cannot be used by more than any one instance of value/quad/root at the same time. Just commonsense really, virtual stack k for model k can only be used for one function evaluation at a time. Obviously there can be at most one instance each of value, root and quad active simultaneously.

r) IFAIL values for D01AJF, D01AEF and C05AZF

Since these iterative techniques may be used inside optimisation or numerical integration procedures, the soft fail option IFAIL = 1 is used. If the Simfit version of these routines is used, a silent exit will occur and failure will not be communicated to users. So it is up to users to be very careful that all calls to quadrature and root finding are consistent and certain to succeed. Default function values of zero are returned on failure.

s) Examples

The test file updown.mod illustrates how to cause a model to swap from one definition to another at critical values of the independent variable, updownup.mod illustrates how to construct a model with two swap-over points, while the series of test files usermodx.tf? give examples of various techniques using sub-models for evaluation, root finding and adaptive quadrature.

t) Synonyms and abbreviations

The following sets of commands are equivalent:-

 sub, minus, subtract
 mul, multiply
 div, divide
 sqrt, squarero, squareroot
 exp, exponent, exponential
 ten, tentothe, tentothepower
 ln, log
 sin, sine
 cos, cosine
 tan, tangent
 asin, arcsin, arcsine
 acos, arccos, arccosine
 atan, arctan, arctangent
 dup, duplicate
 exch, exchange, swap
 del, delete, pop
 abs, absolute
 neg, negative
 min, minimum
 max, maximum
 phi, normal, normalcd, normalcdf
 abserr, epsabs
 relerr, epsrel

3. Special functions

This document describes how to call the special functions commonly used in numerical analysis, statistics, mathematical simulation and data fitting, by one-line commands in Simfit user-supplied models. The options provided at Version 5.4 release 3.039 are described, and a knowledge of the Simfit user-defined model syntax, as described in the files w_readme.f5 and w_readme.f6, is presumed.

Command       NAG    Description
=======       ===    ===========
arctanh(x)    S11AAF Inverse hyperbolic tangent
arcsinh(x)    S11AAF Inverse hyperbolic sine
arccosh(x)    S11AAF Inverse hyperbolic cosine
ai(x)         S17AGF Airy function Ai(x)
dai(x)        S17AJF Derivative of Ai(x)
bi(x)         S17AHF Airy function Bi(x)
dbi(x)        S17AKF Derivative of Bi(x)
besj0(x)      S17AEF Bessel function J0
besj1(x)      S17AFF Bessel function J1
besy0(x)      S17ACF Bessel function Y0
besy1(x)      S17ADF Bessel function Y1
besi0(x)      S18ADF Bessel function I0
besi1(x)      S18AFF Bessel function I1
besk0(x)      S18ACF Bessel function K0
besk1(x)      S18ADF Bessel function K1
phi(x)        S15ABF Normal cdf
phic(x)       S15ACF Normal cdf complement
erf(x)        S15AEF Error function
erfc(x)       S15ADF Error function complement
dawson(x)     S15AFF Dawson integral
ci(x)         S13ACF Cosine integral Ci(x)
si(x)         S13ADF Sine integral Si(x)
e1(x)         S13AAF Exponential integral E1(x)
ei(x)         ...... Exponential integral Ei(x)
rc(x,y)       S21BAF Elliptic integral RC
rf(x,y,z)     S21BBF Elliptic integral RF
rd(x,y,z)     S21BCF Elliptic integral RD
rj(x,y,z,r)   S21BDF Elliptic integral RJ
sn(x,m)       S21CAF Jacobi elliptic function SN
cn(x,m)       S21CAF Jacobi elliptic function CN
dn(x,m)       S21CAF Jacobi elliptic function DN
ln(1+x)       S01BAF ln(1 + x) for x near zero
mchoosen(m,n) ...... Binomial coefficient
gamma(x)      S13AAF Gamma function
lngamma(x)    S14ABF log Gamma function
psi(x)        S14ADF Digamma function, (d/dx)log(Gamma(x))
dpsi(x)       S14ADF Trigamma function, (d^2/dx^2)log(Gamma(x))
igamma(x,a)   S14BAF Incomplete Gamma function
igammac(x,a)  S14BAF Complement of Incomplete Gamma function
fresnelc(x)   S20ADF Fresnel C function
fresnels(x)   S20ACF Fresnel S function
bei(x)        S19ABF Kelvin bei function
ber(x)        S19AAF Kelvin ber function
kei(x)        S19ADF Kelvin kei function
ker(x)        S19ACF Kelvin ker function
cdft(x,m)     G01EBF cdf for t distribution
cdfc(x,m)     G01ECF cdf for chi-square distribution
cdff(x,m,n)   G01EDF cdf for F distribution (m = num, n = denom)
cdfb(x,a,b)   G01EEF cdf for beta distribution
cdfg(x,a,b)   G01EFF cdf for gamma distribution
invn(x)       G01FAF inverse normal
invt(x,m)     G01FBF inverse t
invc(x,m)     G01FCF inverse chi-square
invb(x,a,b)   G01FEF inverse beta
invg(x,a,b)   G01FFF inverse gamma
spence(x)     ...... Spence integral: 0 to x of -(1/y)log|(1-y)|
clausen(x)    ...... Clausen integral: 0 to x of -log(2*sin(t/2))
struveh(x,m)  ...... Struve H function order m (m = 0, 1)
struvel(x,m)  ...... Struve L function order m (m = 0, 1)
kummerm(x,a,b)...... Confluent hypergeometric function M(a,b,x)
kummeru(x,a,b)...... U(a,b,x), b = 1 + n, the logarithmic solution
lpol(x,m,n)   ...... Legendre polynomial of the 1st kind, P_n^m(x),
                     -1 =< x =< 1, 0 =< m =< n
abram(x,m)    ...... Abramovitz function order m (m = 0, 1, 2), x > 0,
                     integral: 0 to infinity of t^m exp( - t^2 - x/t)
debye(x,m)    ...... Debye function of order m (m = 1, 2, 3, 4)
                     (m/x^m)[integral: 0 to x of t^m/(exp(t) - 1)]
fermi(x,a)    ...... Fermi-Dirac integral (1/Gamma(1 + a))[integral:
                     0 to infinity t^a/(1 + exp(t - x))]
heaviside(x,a)...... Heaviside unit impulse function h(x - a)
delta(i,j)    ...... Kronecker delta function
impulse(x,a,b)...... Unit impulse function (small b for Dirac delta)
spike(x,a,b)  ...... Triangular spike unit impulse function
gauss(x,a,b)  ...... Gauss unit impulse function (normal pdf)
sqwave(x,a)   ...... Square wave, amplitude 1, period 2a
rtwave(x,a)   ...... Rectified triangle, amplitude 1, period 2a
mdwave(x,a)   ...... Morse dot wave, amplitude 1, period 2a
stwave(x,a)   ...... Sawtooth wave, amplitude 1, period a
rswave(x,a)   ...... Rectified sine wave, amplitude 1, period pi/a
shwave(x,a)   ...... Rectified half sine, amplitude 1, period 2*pi/a
uiwave(x,a,b) ...... Unit impulse, area 1, period a, width b

Summary

Note that, prior to Version 5.4 release 3.010, only the first eight characters were read of each line, but that has now been increased. Also, to allow users to document their models, from now on all lines starting with a ! within the main model will be ignored and treated as comment lines.

Any of the above commands included as a line in a Simfit model or sub-model simply takes the top stack element as argument and replaces it by the function value. The NAG routines indicated can be consulted for details, as equivalent routines, agreeing very closely with the NAG specifications, are used. The soft fail (IFAIL = 1) options have been used so the simulation will not terminate on error condition, a default value will be returned. Obviously it is up to users to make sure that sensible arguments are supplied, for instance positive degrees of freedom, F or chi-square arguments, etc. To help avoid problems that could arise with arguments outside limits, e.g. probabilities less than zero or greater than one, the command middle (synonym mid) is provided.

Using the command middle to avoid singularities

This command is very useful for avoiding undeflow and overflow, etc. and takes three arguments a, w, b, where a < b. If w =< a then a is placed on the stack, if w > b then b is placed on the stack, otherwise w is placed on the stack. For instance, to obtain an inverse for the normal distribution which requires 0 =< x =< 1, we could use the code

0.001
x
0.999
middle
invn(x)

to avoid nonzero IFAIL returns.

Functions with one argument

The top stack element will be popped and used as an argument, so the routines can be used in several ways. For instance the following code

x
phi(x)
f(1)

would simulate model 1 as a normal cdf, while the code

get(4)
phi(x)
f(3)

would return model three as the normal cdf for whatever was stored in storage location 4.

Functions with two arguments

The top stack element is popped and used as x, while the second is popped and used as a, m, or y, as the case may be. For instance the code

10
x
cdft(x,m)

would place the t distribution cdf with 10 degrees of freedom on the stack, while the code

5
0.975
invt(x,m)

would place the critical value for a two-tailed t test with 5 degrees of freedom at a confidence level of 95% on the stack.

Functions with three or more arguments

The procedure is a simple extension of that described for functions of two arguments. First the stack is prepared as ...u,v,w,z,y,x but, after the function call, it would be ...u,v,w,f(x,y,z). For example, the code

z
y
x
rf(x,y,z)
f(11)

would return model 11 as the elliptic function RF, since f(11) would have been defined as a function of at least three variables. However, the code

get(3)
get(2)
get(1)
rd(x,y,z)
1
add
f(7)

would define model(7) as one plus the elliptic function RD evaluated at whatever was stored in locations 3 (i.e. z), 2 (i.e. y) and 1 (i.e. x).

Test files

Three new test files have been supplied to illustrate these commands:

usermods.tf1: special functions with one argument
usermods.tf2: special functions with two arguments
usermods.tf2: special functions with three arguments

These should be used in program usermod by repeatedly, editing, reading in the edited files, simulating, etc. to explore the options. Users can choose which of the options provided is to be used, by simply uncommenting the desired option and leaving all the others commented. Note that these are all set up for f(1) as a function of one variable and that, by commenting and removing comments so that only one command is active at any one time, the models can be plotted as continuous functions. Alternatively singly calculated values can be compared to tabulated values, which should be indistiguishable if your editing is correct.

4. Operations with vectors

This document describes enhancements to the user-defined model syntax, at version 5.4 release 4.011, designed to increase the scope for model simulation, fitting and plotting where vector arithmetic is involved.

There are techniques to allow initialisation of vectors, i.e. arrays of numerical constants during first pass, commands to accelerate model evaluation, methods to ease the development of models that require vector norms or dot products, and options to simplify repetitive operations such as iterative calculations (i.e. loops).

In particular, 1-line commands to evaluate poynomials or approximate functions by Chebyshev series are described. A knowledge of the Simfit user-defined model syntax described in w_readme.f5, w_readme.f6, and w_readme.f7, is assumed.

The command store(j)

This command is similar to the put(j) command, but there is an important difference; the command put(j) is executed every time the model is evaluated, but the command store(j) is only executed when the model file is parsed for the first time. So store(j) is really equivalent to a data initialisation statement at compile time. For instance, the code

3
store(14)

would initialise store(14) to the value 3. If no further put(14) is used, then storage location 14 would preserve the value 3 for all subsequent calculations in the main model or any sub-model, so that storage location 14 could be regarded as a global constant. Of course any put(14) in the main model or any sub-model would overwrite storage location 14. The main use for the store command is to define special values that are to be used repeatedly for model evaluation, e.g. coefficients for a Chebyshev expansion. For this reason there is another very important difference between put(j) and store(j); store(j) must be preceeded by a literal constant, e.g. 3.2e-6, and cannot be assigned as the end result of a calculation, because storage initialisation is done before calculations.

To summarise: store(j) must be preceded by a numerical value, when it pops this top stack element after copying it into storage location j. So the stack length is decreased by one, to initialise storage location j, but only on the first pass through the model, i.e. when parsing.

The command storef(file)

Since it is tedious to define more than a few storage locations using the command store(j), the command storef(*.*) provides a mechanism for initialising an arbitrary number of storage locations at first pass using contiguous locations. The file specified by the storef command is read and, if it is a Simfit vector file, all the successive components are copied into corresponding storage locations. An example of this is the test model file cheby.mod (and the related data file cheby.dat) which should be run using program usermod to see how a global vector is set up for a Chebshev approximation. Other uses could be when a model involves calculations with a set of fixed observations, such as a time series.

To summarise: the command storef(mydatya) will read the components of any n-dimensional Simfit vector file, mydata, into n successive storage locations starting at position 1, but only when the model file is parsed at first pass. Subsequent use of put(j) or store(j) for j in the range (1,n) will overwrite the previous effect of storef(mydata).

The command poly(x,m,n)

This evaluates m terms of a polynomial by Horner's method of nested multiplication, with terms starting at store(n) and proceeding as far as store(m + n - 1). The polynomial will be of degree m - 1 and it will be evaluated in ascending order. For example, the code

1
store(10)
0
store(11)
-1
store(12)
10
3
x
poly(x,m,n)

will place the value of f(x) = 1 - x^2 on the stack, where x is the local argument. Of course, the contents of the storage locations can also be set by put(j) commands which would then overwrite the previous store(j) command. For instance, the following code

5
put(12)
10
3
2
poly(x,m,n)

used after the previous code, would now place the value 21 on the stack, since f(t) = 1 + 5t^2 = 21, and the argument is now t = 2.

To summarise: poly(x,m,n) evaluates a polynomial of degree m - 1, using successive storage locations n, n + 1, n + 2, ..., n + m - 1, i.e. the constant term is storage location n, and the coefficient of degree m - 1 is storage location m + n - 1. The argument is whatever value is on the top of the stack when poly(x,m,n) is invoked. This command takes three arguments x, m, n off the stack and replaces them by one value, so the stack is decreased by two elements. If there is an error in m or n, e.g. m or n negative, there is no error message, and the value f(x) = 0 is returned.

The command cheby(x,m,n)

The Chebyshev coefficients are first stored in locations n to n + m - 1, then the command cheby(x,m,n) will evaluate a Chebyshev expansion using the Broucke method with m terms. Note that the first term must be twice the constant term since, as usual, only half the constant is used in the expansion. This code, for instance, will return the Chebyshev approximation to exp(x).

 2.532132
 store(20)
 1.130318
 store(21)
 0.271495
 store(22)
 0.044337
 store(23)
 0.005474
 store(24)
 20
 5
 x
 cheby(x,m,n)

Note that, if the numerical values are placed on the stack sequentially, then they obviously must be peeled off in reverse order, as follows.

 2.532132
 1.130318
 0.271495
 0.044337
 0.005474
 store(24)
 store(23)
 store(22)
 store(21)
 store(20)
 20
 5
 x
 cheby(x,m,n)

To summarise: cheby(x,m,n) evaluates a Chebyshev approximation with m terms, using successive storage locations n, n + 1, n + 2, ..., n + m - 1, i.e. twice the constant term is in storage location n, and the coefficient of T(m - 1) is in storage location m + n - 1. The argument is whatever value is on the top of the stack when cheby(x,m,n) is invoked. This command takes three arguments x, m, n off the stack and replaces them by one value, so the stack is decreased by two elements. If there is an error in x, m or n, e.g. x not in (-1,1), or m or n negative, there is no error message, and the value f(x) = 0 is returned. Use test model cheby.mod with program usermod to appreciate this command.

The commands l1norm(m,n), l2norm(m,n) and linorm(m,n)

The Lp norms are calculated for a vector with m terms, starting at storage location n, i.e. l1norm calculates the sum of the absolute values, l2norm calculates the Euclidean norm, while linorm calculates the infinity norm (that is, the largest absolute value in the vector).

It should be emphasised that l2norm(m,n) puts the Euclidean norm on the stack, that is the length of the vector (the square root of the sum of squares of the elements) and not the square of the distance. For example, the code

2
put(5)
-4
put(6)
3
put(7)
4
put(8)
1
put(9)
l1norm(3,5)

would place 9 on the stack, while the command l2norm(5,5) would put 6.78233 on the stack, and the command linorm(5,5) would return 4.

To summarise: these commands take two arguments off the stack and calculate either the sum of the absolute values, the square root of the sum of squares, or the largest absolute value in m successive storage locations starting at location n. The stack length is decreased by one since m and n are popped and replaced by the norm. There are no error messages and, if an error is encountered, a zero value is returned.

The commands sum(m,n) and ssq(m,n)

As there are occasions when it is useful to be able to add up the signed values or the squares of values in storage, these commands are provided. For instance, the code

1
2
3
4
put(103)
put(102)
put(101)
put(100)
100
4
sum(m,n)
f(1)
101
3
ssq(m,n)
f(2)

would assign 10 to function 1 and 29 to function 2.

To summarise: these commands take two arguments off the stack and then replace them with either the sum of m storage locations starting at position n, or the sum of squares of m storage locations starting at position n, so decreasing the stack length by 1.

The command dotprod(l,m,n)

This calculates the scalar product of two vectors of length l which are stored in successive locations starting at positions m and n.

To summarise: The stack length is decreased by 2, as three arguments are consumed, and the top stack element is then set equal to the dot product, unless an error is encountered when zero is returned.

Commands to use mathematical constants

The following commands are provided to faciltate model building.

      pi = 3.141592653589793e+00 ... pi
   piby2 = 1.570796326794897e+00 ... pi divided by two
   piby3 = 1.047197551196598e+00 ... pi divided by three
   piby4 = 7.853981633974483e-01 ... pi divided by four
   twopi = 6.283185307179586e+00 ... pi multiplied by two
 root2pi = 2.506628274631000e+00 ... square root of two pi
 deg2rad = 1.745329251994330e-02 ... degrees to radians
 rad2deg = 5.729577951308232e+01 ... radians to degrees
   root2 = 1.414213562373095e+00 ... square root of two
   root3 = 1.732050807568877e+00 ... square root of three
  eulerg = 5.772156649015329e-01 ... Euler's gamma
lneulerg =-5.495393129816448e-01 ... log (Euler's gamma)

To summarise: these constants are merely added passively to the stack and do not affect any existing stack elements. To use the constants, the necessary further instructions would be required. So, for instance, to transform degrees into radial measure, the code

94.25
deg2rad
multiply

would replace the 94.25 degrees by the equivalent radians.

5. Integers and logical functions

Enhancements to user defined model commands at release 4.016, involving integer and logical functions, and dynamically selected user-defined model functions with arbitrary arguments.

Integer functions

Sometimes integers are needed in models, for instance, as exponents, as summation indices, as logical flags, as limits in do loops, or as pointers in case constructs, etc.

So there are special integer functions that take the top argument off the stack whatever number it is (say x) then replace it by an appropriate integer as follows.

int(x): replace x by the integer part of x
nint(x): replace x by the nearest integer to x
sign(x): replace x by -1 if x < 0, by 0 if x = 0, or by 1 if x > 0

When using integers with simfit models it must be observed that only double precision floating point numbers are stored, and all calculations are done with such numbers, so that 0 actually means 0.0 to machine precision. So, for instance, when using these integer functions with real arguments to create logicals or indices for summation, etc. the numbers on the stack that are to be used as logicals or integers are actually transformed dynamically into integers when required at run time, using the equivalent of nint(x) to generate the appropriate integers. Because of this, you should note that code such as

...
11.3
19.7
1.2
int(x)
2.9
nint(x)
divide

would result in 1.0/3.0 being added to the stack (i.e. 0.3333...) and not 1/3 (i.e 0) as it would for true integer division, leading to the stack

..., 11.3, 19.7, 0.3333333

Logical functions

Logical variables are stored in the global storage vector as either 1.0 (so that nint(x) = 1 = true) or as 0.0 (so that nint(x) = 0 = false). The logical functions either generate logical variables by testing the magnitude of the arbitrary stack value (say x) with respect to zero (to machine precision) or they accept only logical arguments (say m or or n) and return an error message if the stack values are not pre-set to 0.0 or 1.0. Note that logical variables (i.e. Booleans) can be stored using put(i) and retrieved using get(i), so that logical tests of any order of complexity can be constructed.

lt0(x)   replace x by 1 if x < 0, otherwise by 0
le0(x)   replace x by 1 if x =< 0, otherwise by 0
eq0(x)   replace x by 1 if x = 0, otherwise by 0
ge0(x)   replace x by 1 if x >= 0, otherwise by 0
gt0(x)   replace x by 1 if x > 0, otherwise by 0
not(m)   replace m by NOT(m), error if m not 0 or 1
and(m,n) replace m and n by AND(m,n), error if m or n not 0 or 1
or(m,n)  replace m and n by OR(m,n), error if m or n not 0 or 1
xor(m,n) replace m and n by XOR(m,n), error if m or n not 0 or 1

Conditional execution

Using these integer and logical functions in an appropriate sequence interspersed by put(.) and get(.) commands, any storage location (say j) can be set up to test whether any logical condition is true or false. So, the commands if(.) and ifnot(.) are provided to select model features depending on logical variables. The idea is to calculate the logical variables using the integer and logical functions, then load them into storage using put(.) commands. The if(.) and ifnot(.) commands inspect the designated storage locations and return 1 if the storage location has the value 1.0 (to machine precision), or 0 otherwise, even if the location is not 0.0 (to machine precision). The logical values returned model code is executed whereas, if a 0 is returned, the next line is missed out.

if(j)     execute the next line only if storage(j) = 1.0
ifnot(j)  execute the next line only if storage(j) = 0.0

Note that very extensive logical tests and blocks of code for conditional executions, do loops, while and case constructs can be generated by using these logical functions sequentially but, because not all the lines of code will be active, the parsing routines will indicate the number of if(.) and ifnot(.) commands and the resulting potentially unused lines of code. This information is not important for correctly formatted models, but it can be used to check or debug code if required.

Consult the test file if.mod to see how to use logical functions.

Arbitrary functions with arbitrary arguments

The submodels described so far for evaluation, integration, root finding, supplied by the simfit calls to the model evaluation subroutines. However, sometimes it is useful to be able to evaluate a submodel with arbitrary arguments added to the stack, or arguments that are functions of the main arguments. Again, it is useful to be able to evaluate an arbitrary function chosen dynamically from a set of submodels indexed by an integer parameter calculated at run time, rather than read in at compile time when the model is first parsed. So, to extend the user-defined model syntax, the command user1(x,m) is provided. The way this works involves three steps:

1) an integer (m) is put on the stack to denote the required model,
2) calculations are performed to put the argument (x) on the stack, then
3) the user defined model is called and the result placed on the stack.

For instance the code

...
14.7
3
11.3
user1(x,m)

would result in

..., 14.7, 12.5

if the value of submodel number 3 is 12.5 at an argument of 11.3.

Similar syntax is used for functions of two and three variables, i.e.

user1(x,m)
user2(x,y,m)
user3(x,y,z,m)

Clearly the integer m can be literal, calculated or retrieved from storage, but it must correspond to a submodel that has been defined in the sequence of submodels, and the calculated arbitrary arguments x, y, z must be sensible values. For instance the commands

2
x
user1(x,m)

and

value(2)

are equivalent. However the first form is more versatile, as the model number (m, 2 in this case) and argument (x, the dummy value in this case) can be altered dynamically as the result of stack calculations, while the second form will always invoke the case with m = 2 and x = the subroutine argument to the model.

The model file user1.mod illustrates how to use the user1(x,m) command.

6. Using standard expressions

Standard expressions using mathematical notation can be used at any point inside the model definition section of the user defined model file. For example the code

      begin{expression}
      f(1) = [p(1) + p(2)x]/[1.0 + p(3)x + p(4)x^2]
      end{expression}

would define model number 1 as a 2:2 rational function such as are used in steady state enzyme kinetics.
Another example could be the lines

      %
      begin{expression}
      f(1) = p(1)y(1) - p(2)y(1)y(2)
      f(2) = -p(3)y(2) + p(4)y(1)y(2)
      end{expression}
      %
      begin{expression}
      j(1) = p(1) - p(2)y(2)
      j(2) = p(4)y(2)
      j(3) = -p(2)y(1)
      j(4) = -p(3) + p(4)y(1)
      end{expression}
      %

defining the Lotka-Volterra predator-prey differential equations with Jacobian.

7. Summary of commands

The syntax is defined by Bardsley and Prasad (Computers Chem. 21, 71-82, 1997) and in the Simfit reference manual and the Simfit readme files. Program usermod can test codes and check models for errors, and there are example test model files that can be run to illustrate the various techniques.

Commands described in w_readme.f5 (The basic commands)

Command        Stack effect
-------        ------------
x              stack -> stack, x
y              stack -> stack, y
z              stack -> stack, z
add            stack, a, b -> stack, (a + b)
subtract       stack, a, b -> stack, (a - b)
multiply       stack, a, b -> stack, (a*b)
divide         stack, a, b -> stack, (a/b)
p(i)           stack -> stack, p(i)         ...  i can be 1, 2, 3, etc
f(i)           stack, a -> stack  ...evaluate model since now f(i) = a
power          stack, a, b -> stack, (a^b)
squareroot     stack, a -> stack, sqrt(a)
exponential    stack, a -> stack, exp(a)
tentothepower  stack, a -> stack, 10^a
ln (or log)    stack, a -> stack, ln(a)
log10          stack, a -> stack, log(a) (to base ten)
pi             stack -> stack, 3.1415927
sine           stack, a -> stack, sin(a)       ... radians not degrees
cosine         stack, a -> stack, cos(a)       ... radians not degrees
tangent        stack, a -> stack, tan(a)       ... radians not degrees
arcsine        stack, a -> stack, arcsin(a)    ... radians not degrees
arccosine      stack, a -> stack, arccos(a)    ... radians not degrees
arctangent     stack, a -> stack, arctan(a)    ... radians not degrees
sinh           stack, a -> stack, sinh(a)
cosh           stack, a -> stack, cosh(a)
tanh           stack, a -> stack, tanh(a)
exchange       stack, a, b -> stack, b, a
duplicate      stack, a -> stack, a, a
pop            stack, a, b -> stack, a
absolutevalue  stack, a -> stack, abs(a)
negative       stack, a -> stack , -a
minimum        stack, a, b -> stack, min(a,b)
maximum        stack, a, b -> stack, max(a,b)
gammafunction  stack, a -> stack, gamma(a)
lgamma         stack, a -> stack, ln(gamma(a))
normalcdf      stack, a -> stack, phi(a)  integral from -infinity to a
erfc           stack, a -> stack, erfc(a)
y(i)           stack    -> stack, y(i)                Only diff. eqns.
j(i)           stack, a -> stack  J(i-(i/n), (i/n)+1) Only diff. eqns.
***            stack -> stack, ***          ...  *** can be any number

Commands described in w_readme.f6 (global storage and sub-models)

Command           Effect
-------           ------
begin{model(i)}   Start of definition of sub-model i
end{model(i)}     End of definition of sub-model i
put(i)            Cut top stack element into storage location i
get(i)            Copy storage location i to top of stack
get3(i,j,k)       Get one of (i,j,k) following use of command order
epsabs            Re-define absolute error tolerance
epsrel            Re-define relative error tolerance
blim(i)           Re-define bottom limit i
tlim(i)           Re-define top limit i
putpar            Copy parameters into global storage for sub-models
value(i)          Evaluate model i
quad(i)           Integrate model i
convolute(i,j)    Convolute models i and j
root(i)           Estimate a root of model i
value3(i,j,k)     Evaluate one of (i,j,k) following use of command order
order             -1,0,1 depending on values (use before get3 and value3)
middle            Reflect a value back between limits if necessary

Commands described in w_readme.f7 (Special functions)

Command       NAG    Description
-------       ----   -----------
arctanh(x)    S11AAF Inverse hyperbolic tangent
arcsinh(x)    S11AAF Inverse hyperbolic sine
arccosh(x)    S11AAF Inverse hyperbolic cosine
ai(x)         S17AGF Airy function Ai(x)
dai(x)        S17AJF Derivative of Ai(x)
bi(x)         S17AHF Airy function Bi(x)
dbi(x)        S17AKF Derivative of Bi(x)
besj0(x)      S17AEF Bessel function J0
besj1(x)      S17AFF Bessel function J1
besy0(x)      S17ACF Bessel function Y0
besy1(x)      S17ADF Bessel function Y1
besi0(x)      S18ADF Bessel function I0
besi1(x)      S18AFF Bessel function I1
besk0(x)      S18ACF Bessel function K0
besk1(x)      S18ADF Bessel function K1
phi(x)        S15ABF Normal cdf
phic(x)       S15ACF Normal cdf complement
erf(x)        S15AEF Error function
erfc(x)       S15ADF Error function complement
dawson(x)     S15AFF Dawson integral
ci(x)         S13ACF Cosine integral Ci(x)
si(x)         S13ADF Sine integral Si(x)
e1(x)         S13AAF Exponential integral E1(x)
ei(x)         ...... Exponential integral Ei(x)
rc(x,y)       S21BAF Elliptic integral RC
rf(x,y,z)     S21BBF Elliptic integral RF
rd(x,y,z)     S21BCF Elliptic integral RD
rj(x,y,z,r)   S21BDF Elliptic integral RJ
sn(x,m)       S21CAF Jacobi elliptic function SN
cn(x,m)       S21CAF Jacobi elliptic function CN
dn(x,m)       S21CAF Jacobi elliptic function DN
ln(1+x)       S01BAF ln(1 + x) for x near zero
mchoosen(m,n) ...... Binomial coefficient
gamma(x)      S13AAF Gamma function
lngamma(x)    S14ABF log Gamma function
psi(x)        S14ADF Digamma function, (d/dx)log(Gamma(x))
dpsi(x)       S14ADF Trigamma function, (d^2/dx^2)log(Gamma(x))
igamma(x,a)   S14BAF Incomplete Gamma function
igammac(x,a)  S14BAF Complement of Incomplete Gamma function
fresnelc(x)   S20ADF Fresnel C function
fresnels(x)   S20ACF Fresnel S function
bei(x)        S19ABF Kelvin bei function
ber(x)        S19AAF Kelvin ber function
kei(x)        S19ADF Kelvin kei function
ker(x)        S19ACF Kelvin ker function
cdft(x,m)     G01EBF cdf for t distribution
cdfc(x,m)     G01ECF cdf for chi-square distribution
cdff(x,m,n)   G01EDF cdf for F distribution (m = num, n = denom)
cdfb(x,a,b)   G01EEF cdf for beta distribution
cdfg(x,a,b)   G01EFF cdf for gamma distribution
invn(x)       G01FAF inverse normal
invt(x,m)     G01FBF inverse t
invc(x,m)     G01FCF inverse chi-square
invb(x,a,b)   G01FEF inverse beta
invg(x,a,b)   G01FFF inverse gamma
spence(x)     ...... Spence integral: 0 to x of -(1/y)log|(1-y)|
clausen(x)    ...... Clausen integral: 0 to x of -log(2*sin(t/2))
struveh(x,m)  ...... Struve H function order m (m = 0, 1)
struvel(x,m)  ...... Struve L function order m (m = 0, 1)
kummerm(x,a,b)...... Confluent hypergeometric function M(a,b,x)
kummeru(x,a,b)...... U(a,b,x), b = 1 + n, the logarithmic solution
lpol(x,m,n)   ...... Legendre polynomial of the 1st kind, P_n^m(x),
                     -1 =< x =< 1, 0 =< m =< n
abram(x,m)    ...... Abramovitz function order m (m = 0, 1, 2), x > 0,
                     integral: 0 to infinity of t^m exp( - t^2 - x/t)
debye(x,m)    ...... Debye function of order m (m = 1, 2, 3, 4)
                     (m/x^m)[integral: 0 to x of t^m/(exp(t) - 1)]
fermi(x,a)    ...... Fermi-Dirac integral (1/Gamma(1 + a))[integral:
                     0 to infinity t^a/(1 + exp(t - x))]
heaviside(x,a)...... Heaviside unit impulse function h(x - a)
delta(m,n)    ...... Kronecker delta impulse function
impulse(x,a,b)...... Unit impulse function (small b for Dirac delta)
spike(x,a,b)  ...... Unit triangular impulse function
gauss(x,a,b)  ...... Unit Gauss pdf impulse function
sqwave(x,a)   ...... Square wave, amplitude 1, period 2a
rtwave(x,a)   ...... Rectified triangle, amplitude 1, period 2a
mdwave(x,a)   ...... Morse dot wave, amplitude 1, period 2a
stwave(x,a)   ...... Sawtooth wave, amplitude 1, period a
rswave(x,a)   ...... Rectified sine wave, amplitude 1, period pi/a
shwave(x,a)   ...... Sine half-wave, amplitude 1, period 2*pi/2
uiwave(x,a,b) ...... Unit impulse, area 1, period a, width b

Commands described in w_readme.f8 (Operations with vectors)

Command         Effect
-------         ------
store(j)        Initialise global store(j)
storef(file)    Initialise global store(1) to global store(n) from a file
poly(x,m,n)     Evaluate a polynomial
cheby(x,m,n)    Evaluate a Chebyshev expansion
l1norm(m,n)     Evaluate a L_1 norm
l2norm(m,n)     Evaluate a L_2 norm
linorm(m,m)     Evaluate a L_infinity norm
sum(m,n)        Evaluate the sum of vector components
ssq(m,n)        Evaluate the sum of squares of vector components
dotprod(l,m,n)  Evaluate the scalar product of two vectors

Commands to use mathematical constants

 Command   Value                     Comment
 -------   -----                     -------
      pi   3.141592653589793e+00     pi
   piby2   1.570796326794897e+00     pi divided by two
   piby3   1.047197551196598e+00     pi divided by three
   piby4   7.853981633974483e-01     pi divided by four
   twopi   6.283185307179586e+00     pi multiplied by two
 root2pi   2.506628274631000e+00     square root of two pi
 deg2rad   1.745329251994330e-02     degrees to radians
 rad2deg   5.729577951308232e+01     radians to degrees
   root2   1.414213562373095e+00     square root of two
   root3   1.732050807568877e+00     square root of three
  eulerg   5.772156649015329e-01     Euler's gamma
lneulerg  -5.495393129816448e-01     log (Euler's gamma)

Commands described in w_readme.f9 (integers and arbitrary functions)

Note: x, y, z, m and n are the stack values but j is a literal constant, while logicals 1 = true and 0 = false are stored as floating point numbers, to be converted into integers using nint(.) internally when required.

Command        Effect
-------        ------
int(x)         x replaced by the integer part of x
nint(x)        x replaced by the nearest integer to x
sign(x)        x replaced by -1 if x < 0, by 0 if x = 0 or by 1 if x > 0
lt0(x)         x replaced by 1 if x < 0, otherwise by 0
le0(x)         x replaced by 1 if x =< 0, otherwise by 0
eq0(x)         x replaced by 1 if x = 0, otherwise by 0
ge0(x)         x replaced by 1 if x >= 0, otherwise by 0
gt0(x)         x replaced by 1 if x > 0, otherwise by 0
not(m)         m replaced by not(m), error if m not 0 or 1
and(m,n)       m and n replaced by and(m,n), error if m and n not 0 or 1
or(m,n)        m and n replaced by or(m,n), error if m and n not 0 or 1
xor(m,n)       m and n replaced by xor(m,n), error if m and n not 0 or 1
if(j)          execute the next line iff storage(j) = 1.0
ifnot(j)       execute the next line iff storage(j) = 0.0
user1(x,m)     model m evaluated at x
user2(x,y,m)   model m evaluated at x, y
user3(x,y,z,m) model m evaluated at x, y, z

8. Examples

To appreciate the SIMFIT conventions, examine the following test files:-

      deqmod?.tf? ...systems of differential eqns. (require deqpar?.tf?)
      usermod1.tf?...functions of 1 independent variable
      usermod2.tf?...functions of 2 independent variables
      usermod3.tf?...functions of 3 independent variables
      usermodd.tf?...one differential eqn.

The files usermod1.tf? are functions of 1 variable, usermod2.tf? are functions of two variables, usermod3.tf? are functions of 3 variables, while usermodd.tf? are differential equations. Models deqmod?.tf? are models for several functions at the same time, e.g. sets of differential equations as used by deqsol or sqpfit.

For differential equations, the Jacobian must be supplied if you are using program qnfit, but program deqsol will calculate the Jacobian by numerical approximation if you end the model with % then follow with a line starting with a %, instead of the definitions for j(i).

Program usermod can be used to develop, then test a model before trying to use it. The idea is first to examine the test files supplied, then read them into usermod for practise before finally developing your own models. Any file you supply to deqsol, makdat, qnfit, etc. must be formatted exactly like the test files provided.

Model 1 to 6 are first displayed using expressions and then using the less familar reverse Polish notation.

Example 1: straight line

Model as an expression

      %
      Example: p(1) + p(2)*x
      %
      1 equation
      1 variable
      2 parameters
      %
      begin{expression}
      f(1) = p(1) + p(2)x
      end{expression}
      %

Same model in reverse Polish with comments to explain the rules

 %                     start of text defining model indicated by %
Example: user supplied function of 1 variable ... a straight line
.............
p(1) + p(2)*x
.............
 %                     end of text, start of parameters indicated by %
 1 equation            number of equations to define the model
 1 variable            number of variables (or differential equation)
 2 parameters          number of parameters in the model
 %                     end of parameters, start of model indicated by %
 p(1)                  put p(1) on the stack: stack = p(1)
 p(2)                  put p(2) on the stack: stack = p(1), p(2)
 x                     put an x on the stack: stack = p(1), p(2), x
 multiply              multiply top elements: stack = p(1), p(2)*x
 add                   add the  top elements: stack = p(1) + p(2)*x
 f(1)                  evaluate the model f(1) = p(1) + p(2)*x
 %                     end of the model definitions indicated by %

Example 2: single exponential

Model as an expression

      %
      Example: f(x) = p(1)*exp[p(2)*x]
      %
      1 equation
      1 variable
      2 parameters
      %
      begin{expression}
      f(1) = p(1)exp[p(2)x]
      end{expression}
      %

Same model in reverse Polish with comments to explain the rules

 %
Example: user supplied function of 1 variable ... single exponential
                           ....
f(x) = p(1)*exp[p(2)*x]    ....
                           ....
 %
 1 equation            number of equations
 1 variable            number of variables (or differential equation)
 2 parameters          number of parameters in this model
 %
 p(2)                  put p(2) on the stack: stack=p(2)
 x                     put an x on the stack: stack=p(2), x
 multiply              multiply top elements: stack=p(2)*x
 exponential           exponential operator : stack=exp[p(2)*x]
 p(1)                  put p(1) on the stack: stack=exp[p(2)*x], p(1)
 multiply              multiply top elements: stack=p(1)*exp[p(2)*x]
 f(1)                  f(1) = p(1)*exp[p(2)*x]
 %

Example 3: normal integral

Model as an expression

      %
      Example: normal integral
      %
      1 equation
      1 variable
      3 parameters
      %
      begin{expression}
      f(1) = p(3)normalcdf[(x - p(1))/p(2)]
      end{expression}
      %

Same model in reverse Polish with comments to explain the rules

 %
Example: user supplied function of 1 variable ... normal integral
         i.e.,  integral from -infinity to x of
p(3)/[p(2)*sqrt(2*pi)]*exp{(-1/2)*[(u - p(1))/p(2)]^2}
  %
 1 equation         number of equations
 1 variable         number of variables (or differential equation)
 3 parameters       number of parameters in this model
 %
 x                  stack=x
 p(1)               stack=x, p(1)
 subtract           stack=x - p(1)
 p(2)               stack=x - p(1), p(2)
 divide             stack=[x - p(1)]/p(2)
 normalcdf          stack=integral
 p(3)               stack=integral, p(3)
 multiply           stack=p(3)*integral
 f(1)               evaluate model
 %

Example 4: capillary diffusion

Model as an expression

      %
      Example: f(x) = p(1)*erfc[x/(2*sqrt(p(2))]
      %
      1 equation
      1 variable
      2 parameters
      %
      begin{expression}
      f(1) = p(1)erfc[x/(2sqrt(p(2))]
      end{expression}
      %

Same model in reverse Polish

 %
Example: user supplied function of 1 variable ... capillary diffusion
f(x) = p(1)*erfc[x/(2*sqrt(p(2))]
 %
 1 equation
 1 variable
 2 parameters
 %
 x
 p(2)
 squareroot
 2
 multiply
 divide
 erfc
 p(1)
 multiply
 f(1)
 %

Example 5: damped simple harmonic motion

Model as an expression

      %
      Example: f(x) = p(4)*exp[-p(3)*x]*cos[p(1)*x - p(2)]
      %
      1 equation
      1 variable
      4 parameters
      %
      begin{expression}
      f(1) = p(4)exp[-p(3)x]cos[p(1)x - p(2)]
      end{expression}
      %

Same model in reverse Polish.

 %
Example: user supplied function of 1 variable ... damped SHM
         Damped simple harmonic motion in the form
         f(x) = p(4)*exp[-p(3)*x]*cos[p(1)*x - p(2)]
         where p(i) >= 0
 %
 1 equation
 1 variable
 4 parameters
 %
 p(1)
 x
 multiply
 p(2)
 subtract
 cosine
 p(3)
 x
 multiply
 negative
 exponential
 multiply
 p(4)
 multiply
 f(1)
 %

Example 6: Lotka-Voterra predator-prey differential equations

Model using expressions

      %
      Example: Lotka-Volterra predator-prey equations
      differential equations: f(1) = dy(1)/dx
                                   = p(1)*y(1) - p(2)*y(1)*y(2)
                              f(2) = dy(2)/dx
                                   = -p(3)*y(2) + p(4)*y(1)*y(2)
                    jacobian: j(1) = df(1)/dy(1)
                                   = p(1) - p(2)*y(2)
                              j(2) = df(2)/dy(1)
                                   = p(4)*y(2)
                              j(3) = df(1)/dy(2)
                                   = -p(2)*y(1)
                              j(4) = df(2)/dy(2)
                                   = -p(3) + p(4)*y(1)
          initial condition: y0(1) = p(5), y0(2) = p(6)
      %
      2 equations
      differential equation
      6 parameters
      %
      begin{expression}
      f(1) = p(1)y(1) - p(2)y(1)y(2)
      f(2) = -p(3)y(2) + p(4)y(1)y(2)
      end{expression}
      %
      begin{expression}
      j(1) = p(1) - p(2)y(2)
      j(2) = p(4)y(2)
      j(3) = -p(2)y(1)
      j(4) = -p(3) + p(4)y(1)
      end{expression}
      %

The coding for the model is now finished but, optionally, parameter starting values, curve fitting limits and the range of integration can be appended. If these are not supplied program DEQSOL will request an initialisation file.

 begin{limits}
 0  1.0  3
 0  1.0  3
 0  1.0  3
 0  1.0  3
 0  1.0  3
 0  0.5  3
 end{limits}
 begin{range}
 121
 0
 10
 end{range}

Same model in reverse Polish.

 %
 Example of a user supplied pair of differential equations
 file: deqmod2.tf2 (typical parameter file deqpar2.tf2)
 model: Lotka-Volterra predator-prey equations
 differential equations: f(1) = dy(1)/dx
                              = p(1)*y(1) - p(2)*y(1)*y(2)
                         f(2) = dy(2)/dx
                              = -p(3)*y(2) + p(4)*y(1)*y(2)
              jacobian: j(1) = df(1)/dy(1)
                             = p(1) - p(2)*y(2)
                        j(2) = df(2)/dy(1)
                             = p(4)*y(2)
                        j(3) = df(1)/dy(2)
                             = -p(2)*y(1)
                        j(4) = df(2)/dy(2)
                             = -p(3) + p(4)*y(1)
    initial condition: y0(1) = p(5), y0(2) = p(6)
 Note: the last parameters must be y0(i) in differential equations
 %
 2 equations
 differential equation
 6 parameters
 %
 y(1)
 y(2)
 multiply
 duplicate
 p(2)
 multiply
 negative
 p(1)
 y(1)
 multiply
 add
 f(1)
 p(4)
 multiply
 p(3)
 y(2)
 multiply
 subtract
 f(2)
 %
 p(1)
 p(2)
 y(2)
 multiply
 subtract
 j(1)
 p(4)
 y(2)
 multiply
 j(2)
 p(2)
 y(1)
 multiply
 negative
 j(3)
 p(4)
 y(1)
 multiply
 p(3)
 subtract
 j(4)
 %

Example 7: convolution integral

Model using expressions

%
convolution integral: from 0 to x of f(u)*g(x - u) du, where
f1(t) = f(t) = exp(-p(1)*t)
f2(t) = g(t) = [p(2)^2]*t*exp(-p(2)*t)
f3(t) = f*g  = f1*f2
--------------------------------------------------------------
This demonstrates how to define 2 equations as sub-models, using
expressions A, B, C  to communicate parameters to the sub-models,
and the command convolute(1,2) to integrate sub-models 1 and 2
(by adaptive quadrature) from blim(1) = 0 to t = tlim(1) = x.
Precision of D01AJF quadrature is controlled by epsabs and epsrel
and blim(1) and tlim(1) must be used for the convolution limits
which, in this case are 0 to x, where x > 0 by assumption.
Note that usually extra parameters must be supplied if it wished
to normalise so that the integral of f or g or f*g is specified
(e.g.,equals 1) over the total range of possible integration.
This must often be done, e.g., if g(.) is a density function.
The gamma distribution normalising factor p(2)**2 is stored as
C in this example to avoid unnecessary re-calculation.
 %
 3 equations
 1 variable
 2 parameters
 %
 begin{expression}
 A = p(1)
 B = p(2)
 C = p(2)*p(2)
 end{expression}
 1
 x
 user1(x,m)
 f(1)
 2
 x
 user1(x,m)
 f(2)
 0.0001
 epsabs
 0.001
 epsrel
 0
 blim(1)
 x
 tlim(1)
 convolute(1,2)
 f(3)
 %
begin{model(1)}
 %
Example: exponential decay, exp(-p(1)*x)
 %
 1 equation
 1 variable
 0 parameter
 %
 begin{expression}
 f(1) = exp(-A*x)
 end{expression}
 %
end{model(1)}
begin{model(2)}
 %
Example: gamma density of order 2
 %
 1 equation
 1 variable
 0 parameters
 %
 begin{expression}
 f(1) = C*x*exp(-B*x)
 end{expression}
 %
end{model(2)}
Advice: consult w_readme.f6 for details of how to define sub models

Same model using reverse Polish

 %
An example using two sub-models for a convolution integral f*g
=================================================================
This demonstrates how to define 2 equations as sub-models, using
the command putpar to communicate parameters to the sub-models,
and the command convolute(1,2) to integrate sub-models 1 and 2
(by adaptive quadrature) from blim(1) = 0 to t = tlim(1) = x.
Precision of D01AJF quadrature is controlled by epsabs and epsrel
and blim(1) and tlim(1) must be used for the convolution limits
which, in this case are 0 to x, where x > 0 by assumption.
.................................................................
integral: from 0 to x of f(u)*g(x - u) du, where
f(t) = exp(-p(1)*t)
g(t) = [p(2)^2]*t*exp(-p(2)*t)
.................................................................
Note that usually extra parameters must be supplied if it wished
to normalise so that the integral of f or g or f*g is specified
(e.g.,equals 1) over the total range of possible integration.
This must often be done, e.g., if g(.) is a density function.
The gamma distribution normalising  factor p(2)**2 is stored in
this example to avoid unnecessary re-calculation.
 %
 1 equation
 1 variable
 2 parameters
 %
 putpar
 p(2)
 p(2)
 multiply
 put(1)
 0.0001
 epsabs
 0.001
 epsrel
 0
 blim(1)
 x
 tlim(1)
 convolute(1,2)
 f(1)
 %
begin{model(1)}
 %
Example: exponential decay, exp(-p(1)*x)
 %
 1 equation
 1 variable
 1 parameter
 %
 p(1)
 x
 multiply
 negative
 exponential
 f(1)
 %
end{model(1)}
begin{model(2)}
 %
Example: gamma density of order 2
 %
 1 equation
 1 variable
 2 parameters
 %
 p(2)
 x
 multiply
 negative
 exponential
 x
 multiply
 get(1)
 multiply
 f(1)
 %
end{model(2)}

Example 8: Rosenbrock's function

Model using expressions

      %
      Example: Rosenbrock's 2-dimensional function for optimisation
      f(1) = 100(y - x^2)^2 + (1 - x)^2
      f(2) = g(1)
           = d(f(1))/dx
           = -400x(y - x^2) - 2(1 - x)
      f(3) = g(2)
           = d(f(1))/dy
           = 200(y - x^2)
      %
      3 equations
      2 variables
      0 parameters
      %
      begin{expression}
      A = y - x^2
      B = 1 - x
      f(1) = 100*A^2 + B^2
      f(2) = -400A*x - 2B
      f(3) = 200A
      end{expression}
      %

Same model in reverse Polish

 %
Rosenbrock's 2-dimensional function for optimisation
f(1) = 100(y - x^2)^2 + (1 - x)^2
f(2) = g(1)
     = d(f(1))/dx
     = -400x(y - x^2) - 2(1 - x)
f(3) = g(2)
     = d(f(1))/dy
     = 200(y - x^2)
%
 3 equations
 2 variables
 0 parameters
 %
 y
 x
 x
 multiply
 subtract
 put(1)
 1.0
 x
 subtract
 put(2)
 get(1)
 get(1)
 multiply
 100.0
 multiply
 get(2)
 get(2)
 multiply
 add
 f(1)
 get(1)
 x
 multiply
 -400.0
 multiply
 get(2)
 -2.0
 multiply
 add
 f(2)
 get(1)
 200.0
 multiply
 f(3)
 %

Example 9: Helix in parameteric form

Model using expressions

      %
      Example: X = A*cos(t), Y = B*sin(t), Z = C*t
               where: t = x, A = p(1), B = p(2), C = p(3)
                      and X(t) = f(1), Y(t) = f(2), Z(t) = f(3)
      %
      3 equations
      1 variable
      3 parameters
      %
      begin{expression}
      f(1) = p(1)cos(x)
      f(2) = p(2)sin(x)
      f(3) = p(3)x
      end{expression}
      %

Same model in reverse Polish

 %
Example: the helix
         X = A*cos(t), Y = B*sin(t), Z = C*t
         where: t = x, A = p(1), B = p(2), C = p(3)
                and X(t) = f(1), Y(t) = f(2), Z(t) = f(3)
 %
 3 equations
 1 variable
 3 parameters
 %
 x
 cos
 p(1)
 multiply
 f(1)
 x
 sin
 p(2)
 multiply
 f(2)
 x
 p(3)
 multiply
 f(3)
 %

Simfit home page