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.

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

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.

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:

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(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.

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.

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

Test files

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).

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

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.

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

5. Integers and logical functions

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

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.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.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.

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.

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.