A model can be solved numerically from the initial state by calling
run(), and the output will be a timeplot, trajectory or
table. Next to the graphics output, run() returns the final
state attained by the simulation (or all data when
table=TRUE). The former can be helpful if one wants to
continue from the previous state (e.g.,
f<-run(); f<-run(state=f)). The full definition of
run() is:
run <- function(tmax=100, tstep=1, state=s, parms=p, odes=model, ymin=0, ymax=NULL, log="", xlab="Time", ylab="Density", tmin=0, draw=lines, times=NULL, show=NULL, arrest=NULL, events=NULL, after=NULL, tweak=NULL, timeplot=TRUE, traject=FALSE, table=FALSE, add=FALSE, legend=TRUE, solution=FALSE, delay=FALSE, lwd=2, col="black", pch=20, ...) {}run() calls the ode() function from the
deSolve package. Additional arguments, like
main and sub, are passed on to
ode() and plot() via the ....
The phase plane function plane() sets up a space with
the first variable on the horizontal axis, and the second on the
vertical axis. The full definition of plane() is:
plane <- function(xmin=-0.001, xmax=1.05, ymin=-0.001, ymax=1.05, xlab="", ylab="", log="", npixels=500, state=s, parms=p, odes=model, x=1, y=2, time=0, grid=5, show=NULL, addone=FALSE, portrait=FALSE, vector=FALSE, add=FALSE, legend=TRUE, zero=TRUE, lwd=2, col="black", pch=20, ...) {}Additional arguments ... are passed on to
run() (for the phase portrait) and to plot()
(e.g., for main and sub). Note that
plane() calls the “vectorized” R-function
outer(), which implies that if one calls functions in the
ODEs they should also be vectorized, e.g., one should use
pmax() instead of max(). There is an
extension, cube.R, for 3-dimensional nullclines.
The function newton() finds a steady state from a nearby
initial state, and can report the Jacobi matrix with its eigenvalues and
eigenvectors. The full definition of newton() is:
newton <- function(state=s, parms=p, odes=model, time=0, positive=FALSE, jacobian=FALSE, vector=FALSE, plot=FALSE, silent=FALSE, addone=FALSE, ...) { }newton() calls the function steady() from
the rootSolve package (which calls stode()).
Additional arguments, ..., are passed on to both of them.
newton() needs an initial state close to an equilibrium
point.
The function continue() continues a steady state by
changing a “bifurcation" parameter defined by the horizontal axis of the
bifurcation diagram. The full definition of continue()
is:
continue <- function(state=s, parms=p, odes=model, step=0.01, x=1, y=2, time=0, xmin=0, xmax=1,ymin=0, ymax=1.05, xlab="", ylab="", log="", col=c("red","black","blue"), lwd= c(2,1,1), addone=FALSE, positive=FALSE, nvar=FALSE, add=FALSE, ...) { }continue() calls the function steady()from
the rootSolve package (additional arguments (...) are
passed on), and needs an initial state close to an equilibrium point.
Note that there is much more proper software for bifurcation analysis
like XPPAUT or MatCont, which reports the type of bifurcations
encountered, and automatically continues all branches of branch points.
Additionally, deBif is a Grind-like R-package allowing for
interactive bifurcation analysis using a shiny-based interface.
The function fit() fits a model to data by non-linear
parameter estimation. The output is an object containing the estimated
parameters, the summed squared residuals, confidence ranges, and
correlations (see the modFit() manual). The data and the
model behavior for its best fit parameters are depicted in a timeplot.
Its full definition is:
fit <- function(datas=data, state=s, parms=p, odes=model, free=NULL, who=NULL, differ=NULL, fixed=NULL, tmin=0, tmax=NULL, ymin=NULL, ymax=NULL, log="", xlab="Time", ylab="Density", bootstrap=0, show=NULL, fun=NULL, costfun=cost, logpar=FALSE, lower=-Inf, upper=Inf, initial=FALSE, add=FALSE, timeplot=TRUE, legend=TRUE, main=NULL, sub=NULL, pchMap=NULL, ...) { }fit() calls the function modFit() from the
FME package (which calls modCost()).
Additional arguments, ..., are passed on to both of them,
and to run() and ode().
Finally the internal function timePlot() can be used to
plot data and is defined as:
timePlot <- function(data, tmin=0, tmax=NULL, ymin=0, ymax=NULL, log="", xlab="Time", ylab="Density", show=NULL, legend=TRUE, draw=lines, lwd=2, add=FALSE, main=NULL, sub=NULL, colMap=NULL, pchMap=NULL, ...) { }These functions have many options (arguments). Fortunately most of them have good default values, and can typically be omitted. Since many options are the same in the five Grind functions, we list them alphabetically, giving their default values, and a short explanation.
add=FALSE defines whether a new plot should be created
when drawing new graphics. When add=TRUE Grind assumes the
definition of the axes have not been changed. Used by
continue(), run(), and
plane().addone=FALSE allows for a “zero" on logarithmic axis
(like the log1p() function in R). Setting
addone=TRUE and xmin or ymin to
one (or 0.99) allows one to include steady states in which one of the
variables is zero in logarithmic nullclines and parameter continuations.
Used by plane(), continue(), and
newton().after=NULL defines commands to be executed after each
time step of a numerical integration, e.g.,
after="parms[1]<-rnorm(1,1,0.1)", draws a new random
value of the first parameter after each time step (as defined by
tstep). Used by run() and
fit().arrest=NULL defines a vector of values, or parameter
names, defining time points where the integrator should stop, and report
the current state (i.e., these time points are added to the
times vector of ode()). This can be helpful
for running (and fitting) models with if statements, or
other discontinuous time points, e.g.,
arrest=c("T1","T2")). Used by run() and
fit().bootstrap=0 defines the number of samples to be taken
randomly from the data (with replacement). This prints a summary and
adds an element bootstrap to the modFit list,
containing a matrix with all parameter estimates. Use
pairs(f$bootstrap) to see the correlations between the
estimates. Used by fit().col="black" in and defines the color of
trajectories.col=c("red","black","blue") in defines the color of
stable, neutral, and unstable steady states, respectively.
Alternatively, when nvar=TRUE one can use this option to
color a steady state by its number of non-zero variables (e.g.,
col=colors would do if you like Grind’s default
colors).colMap=NULL is vector of integers used to re-map the
colors, e.g., colMap=c(3,2,1) reverts the order of the
first three colors. Note that grind.R defines its own color
table (of dark colors that print well). Used by continue(),
plane() and run().costfun=cost allows one to redefine the cost-function
measuring the distance between the model and the data. This can be
useful when different data sets need different models. The default
cost-function loops over the various data sets, and calls
run() for each of them. That call can easily be adapted for
each data set (the index of the loop is called iset). Used
by fit().datas=data in fit() defines the name of
the data frame containing the data, or defines a list of data frames.
Used by fit().delay=FALSE tells Grind whether or not to call the DDE
solver from the deSolve package. The time delay(s) in the
model are to be defined by calling the lagvalue() function
from the deSolve package. This can only be used in
run() and fit() and should also not be used in
combination with phase plane analysis, nor with searching steady states
(newton(), continue()). Used by run() and
fit().differ=NULL defines the names of the free parameters
that differ between the data sets and need to be fitted separately.
differ can also be a named list containing the individual
guesses for each data set. (One can use
makelist(differ,state,parms,nsets) to set up such a list).
Used by fit().draw=lines draws time plots as continuous lines. The
alternative is draw=points. This tells
timePlot() to use the R-functions lines() or
points() to plot the columns of the data frame. Used by
run() and fit().eps was a shortcut in plane() to include
or exclude a zero-axis in the drawing of nullclines. The option has been
removed because the defaults of xmin and ymin
in plane() were set to -0.001, which means
that a nullcline on the axes will be drawn.events is an option that can be passed on to
deSolve() to handle events during numerical integration.
Type events? for help. Used by run().fixed=NULL defines a named list of the non-free
parameters that differ between the data sets and have known fixed
values. (One can use makelist(fixed,state,parms,nsets) to
set up such a list). Used by fit().free=NULL defines the names of the ‘free’ parameters to
be fitted. By default free equals
c(names(state), names(parms)).fit().fun=NULL defines a function to transform the data and
(numerical) solution before fitting. Used by fit().grid=5 defines the number of grid points for which the
vector field or phase portrait is drawn. Used by
plane().initial=FALSE allows one to read the initial condition
from the data (and not estimate it). Still used by fit()
but will disappear.jacobian=FALSE defines whether or not
newton() should print the Jacobian. Used by
newton().log="" can be used to define logarithmic axis (this is
passed on to the R-function plot()). log="xy"
makes both axes logarithmic. Used by continue(),
plane() and run().logpar=FALSE can be used to fit the logarithm of the
parameters, rather than their true values. This can be very useful when
the parameters have very different values. Grind will report the true
values that are estimated, but modFit() returns the
logarithmic values. Used by fit().lower=-Inf defines the lower limit of the free
parameters. One can set lower=0 to keep all free parameters
positive. Otherwise, the order of values in lower should be the same as
the names of the parameters in free. This becomes
complicated when there are shared and different free parameters (as
defined by free and differ). The order in
lower should then first be the parameters in
free and then those in differ. The latter are
automatically repeated for every additional data set if the length of
lower equals the sum of the lengths of free
and differ. Thus, the length of lower is
either 1, length(free)+length(differ) or
length(free)+n*length(differ), where n is the
number of data sets. Used by fit().lwd=2 sets the line-width of graphs. Used by
continue(), plane() and
run().main=NULL allows one to put a title at the top of the
graph (this is passed on to the R-function plot() by all
Grind functions). In fit() main can be a vector for the
titles of all plots when multiple data sets are being fitted. Note that
titles are set in a plain font (to set that back to R’s default bold
face, change the two font lines in the grind.R
file). Used by continue(), fit(),
plane() and run().npixels=500 defines the resolution of the phase space
in plane(). Setting it to a lower value speeds up the
drawing of nullclines. Used by plane().nvar=FALSE asks continue() to color steady
state by their stability. By default stable steady state are depicted as
heavy red lines and unstable steady states by thin blue lines (which can
be changed with the col and lwd options of ).
If nvar=TRUE the states are colored by their number of
non-zero variables (use col to define these colors), and
their stability remains indicated by the line width. Used by
continue().odes=model defines the name the model. Used by all
Grind’s functions.parms=p define the name of the parameter vector. Used
by all Grind’s functions.pch=20 defines the symbol to be plotted at the start of
trajectories. Used by plane() and run().pchMap=NULL can be used to re-map the symbols, e.g.,
pchMap=c(3,2,1) reverts the order of the first three
R-symbols (see pch in the R-function
points()). Used by . Used by run().plot=FALSE defines whether or not newton()
depict the steady state by a symbol in the phase plane. Used by
newton().portrait=FALSE defines whether or not
plane() should include a phase portrait. Used by
plane().positive=FALSE, setting positive=TRUE
restricts the search of newton() and
continue() to positive steady states only. Used by
newton() and continue().silent=FALSE. By default newton() returns
a steady state, while printing the eigenvalues. After setting
silent=TRUE, newton() does not print anything,
and returns a list containing the state, the Jacobian, the eigenvalues
and the eigenvectors. Use
max(Re(newton(s,silent=TRUE)$values)) to retrieve the
largest eigenvalue. Used by newton().show=NULL defines the variables appearing in a time
plot, fit, or phase plane. By default all variables are shown. By
explicitly providing a list of names one can define subsets to be
plotted, e.g., show=c("P","Q"). Used by run()
and plane().solution=FALSE defines whether or not the model
provides time derivatives (default), or a full solution. This can only
be used in run() and fit() and should
obviously not be used in combination with phase plane analysis, nor with
searching steady states. Models returning a solution obey the same
format as the ODE models required by deSolve (but can
return a vector of values and not necessarily a list). Used by
run() and fit().state=s define the name of the state vector. Used by
all Grind’s functions.step=0.01 defines the maximum change of the bifurcation
parameter in a bifurcation diagram. When the axis is linear the
parameter is increased, or decreased, in steps not exceeding
step \(\times\)
xmax. When the axis is logarithmic the parameters is
maximally multiplied by \(1+\)step.
continue() will decrease the step size to maximally
step/100 when it looses the steady state. Used by
continue().sub=NULL allows one to put a subtitle at the bottom of
the graph (this is passed on to the R-function plot()).
Note that titles are set in a plain font (to set that back to bold face,
change the two font lines in the grind.R
file). Used by continue(), plane() and
run().table=FALSE asks run() to return a table
with the values reported by the integrator. Used by
run().time=0 defines the time point for which nullclines are
computed and steady states are computed (for non-autonomous ODEs). Used
by continue(), newton() and
plane().times=NULL can be used to define the time points at
which the numerical integrator generates output (even though this is
typically done by setting tstep). Used by
run().timeplot=TRUE asks run() to make a time
plot. Used by run() and fit().tmax=100 set the integration time. Used by
run() and fit().tmin=0 allows one to start a specific time point (which
can be convenient when a run is continued). Used by
run().traject=FALSE asks run() to draw
trajectory in an existing phase plane. Used by run().tstep=1 set the reporting interval of numerical
integrations. One can also provide a vector of time points where the
intergrator should provide output with the times option
(see the ode() manual). Used by run() and
fit().tweak=NULL allows one to modify the data delivered by
run(). For instance one can add columns that can be fitted
to data, e.g., tweak<-"nsol$T<-nsol$R+nsol$N" to sum
\(R\) and \(N\) into a new column \(T\). One can also transform predictions
from the model before they are fitted to data that is already
transformed. Used by run() and fit().upper=Inf defines the upper limit of the free
parameters. One can set upper=10 to give all free
parameters a maximum value of 10. Otherwise, the order of values in
upper should be the same as the names of the parameters in
free. This becomes complicated when there are shared and
different free parameters (as defined by free and
differ). The order in upper should then first
be the parameters in free and then those in
differ. The latter are automatically repeated for every
additional data set if the length of upper equals the sum
of the lengths of free and differ. Thus, the
length of upper is either 1,
length(free)+length(differ) or
length(free)+n*length(differ), where n is the
number of data sets. Used by fit().x=1 define the variable on the horizontal axis of phase
planes and bifurcation diagrams. One can also use the names of the
variables to define the axes, e.g., x="R". Used by
continue() and plane().xlab="Time" or xlab=variable allows one to
redefine the label of a horizontal axes. Used by
continue(), plane() and
run().xmax=1.1 define the maximum of vertical axes of time
plots phase planes and bifurcation diagrams. Used by
continue() and plane().xmin define the origin of the horizontal axis of phase
planes and bifurcation diagrams. Used by continue() where
xmin=0 by default, and plane() where
xmin=-0.001 is the default.y=2 define the variable on the vertical axis of phase
planes and bifurcation diagrams. One can also use the names of the
variables to define the axes, e.g., y="N". Used by
continue() and plane().ylab="Density" or ylab=variable allows one
to redefine the label of a vertical axes. Used by
continue(), plane() and
run().ymax=1.1 define the maximum of the vertical axis of
time plots, phase planes and bifurcation diagrams. Used by . Used by
continue(), plane() and
run().ymin define the origin of the vertical axis of phase
planes and bifurcation diagrams. Used by continue() where
ymin=0 by default, and plane() where
ymin=-0.001 is the default.vector=FALSE defines whether or not
plane() should include a vector field. Used by
plane().vectors=FALSE defines whether or not
newton() should print the eigenvectors. Used by
newton().zero=TRUE draws the phase plane for all variables other
than x and y set to zero (important when
drawing nullclines of variables not appearing on the axes). Used by
plane().... can be used to define parameters that are passed on
to other functionsCheck Grind’s homepage for more info.