Solve a replicate set of difference (or recurrence) equation by
iterating it a number of times. This is a wrapper around
difeq that does not (yet) do anything clever to
avoid many allocations.
difeq_replicate(n, y, steps, target, parms, ..., n_out = 0L,
n_history = 0L, grow_history = FALSE, return_history = n_history > 0,
dllname = "", parms_are_real = TRUE, ynames = NULL, outnames = NULL,
return_by_column = TRUE, return_initial = TRUE, return_step = TRUE,
return_output_with_y = TRUE, restartable = FALSE,
return_minimal = FALSE, as_array = TRUE)Number of replicates. It is an error to request zero replicates.
The initial state of the system. Must be either a
numeric vector or a list of numeric vectors. If the
latter, it must have length n.
A vector of steps to return the system at. The
first step is taken as step zero, and the solution will
be recorded at every other step in the vector. So to step a
system from time zero to times 1, 2, 3, ..., n use 0:n. Must be
integer values and will be passed through as.integer
(which may truncate or otherwise butcher non-integer values).
The target function to advance. This can either be
an R function taking arguments n, i, t, y, parms or be a
scalar character with the name of a compiled function with
arguments size_t n, size_t step, double time, const double
*y, double *dydt, size_t n_out, double *output void *data.
Parameters to pass through to the difference function
Dummy arguments - nothing is allowed here, but this means that all further arguments must be specified by name (not order) so I can easily reorder them later on.
The number of output variables (in addition to the
difference equation variables). If given, then an R function
must return an attribute output with the output
variables.
The number of iterations of history to save during the simulation. By default, no history is saved.
Logical indicating if history should be grown
during the simulation. If FALSE (the default) then when
history is used it is overwritten as needed (so only the most
recent n_history elements are saved. This may require
some tuning so that you have enough history to run your
simulation (i.e. to the longest delay) or an error will be
thrown when it underflows. The required history length will
vary with your delay sizes and with the timestep for dopri. If
TRUE, then history will grow as the buffer is exhausted.
The growth is geometric, so every time it reaches the end of the
buffer it will increase by a factor of about 1.6 (see the
ring documentation). This may consume more memory than
necessary, but may be useful where you don't want to care about
picking the history length carefully.
Logical indicating if history is to be
returned. By default, history is returned if n_history
is nonzero.
Name of the shared library (without extension) to
find the function func in the case where func
refers to compiled function.
Logical, indicating if parms should
be treated as vector of doubles by func (when it is a
compiled function). If TRUE (the default), then
REAL(parms), which is double* is passed through.
If FALSE then if params is an externalptr type
(EXTPTRSXP) we pass through the result of
R_ExternalPtrAddr, otherwise we pass params
through unmodified as a SEXP. In the last case, in your
target function you will need to include <Rinternals.h>,
cast to SEXP and then pull it apart using the R API (or
Rcpp).
Logical, indicating if the output should be named
following the names of the input vector y.
Alternatively, if ynames is a character vector of the
same length as y, these will be used as the output names.
An optional character vector, used when
n_out is greater than 0, to name the model output matrix.
Logical, indicating if the output should be
returned organised by column (rather than row). This incurs a
slight cost for transposing the matrices. If you can work with
matrices that are transposed relative to deSolve, then
set this to FALSE.
Logical, indicating if the output should
include the initial conditions. Specifying FALSE avoids
binding this onto the output.
Logical, indicating if a row (or column if
return_by_column is TRUE) representing step is included.
Logical, indicating if the output
should be bound together with the returned matrix y (as
it is with deSolve). If FALSE, then output will
be returned as the attribute output.
Logical, indicating if the problem should be
restartable. If TRUE, then the return value of a
simulation can be passed to difeq_restart to continue the
simulation after arbitrary changes to the state or the
parameters. Note that this is really only useful for delay
difference equations where you want to keep the history but make
changes to the parameters or to the state vector while keeping
the history of the problem so far.
Shorthand option - if set to TRUE
then it sets all of return_by_column,
return_initial, return_time,
return_output_with_y to FALSE
(Defunct) Logical, indicating if the output should
be converted into an array. If TRUE then res[, ,
i] will contain the i'th replicate, if FALSE then
res[[i]] does instead. If both as_array and
restartable are TRUE, then the attributes
ptr and restart_data will be present as a
list of restarting information for difeq_continue,
though using these is not yet supported.
It is not currently possible to replicate over a set of parameters at once yet; the same parameter set will be used for all replications.
The details of how replication is done here are all considered
implementation details and are up for change in the future - in
particular if the models are run in turn or simultaneously (and
the effect that has on the random number stream). Logic around
naming output may change in future too; note that varying names in
the y here will have some unexpected behaviours.
# Here is a really simple equation that does a random walk with
# steps that are normally distributed:
rhs <- function(i, y, p) y + runif(1)
y0 <- 1
t <- 0:10
p <- 5
dde::difeq_replicate(10, y0, t, rhs, p)
#> , , 1
#>
#> [,1] [,2]
#> [1,] 0 1.000000
#> [2,] 1 1.080750
#> [3,] 2 1.915083
#> [4,] 3 2.515844
#> [5,] 4 2.673053
#> [6,] 5 2.680452
#> [7,] 6 3.146845
#> [8,] 7 3.644623
#> [9,] 8 3.934390
#> [10,] 9 4.667272
#> [11,] 10 5.439794
#>
#> , , 2
#>
#> [,1] [,2]
#> [1,] 0 1.000000
#> [2,] 1 1.874601
#> [3,] 2 2.049541
#> [4,] 3 2.083783
#> [5,] 4 2.404168
#> [6,] 5 2.806497
#> [7,] 6 3.002166
#> [8,] 7 3.405705
#> [9,] 8 3.469366
#> [10,] 9 3.858067
#> [11,] 10 4.833615
#>
#> , , 3
#>
#> [,1] [,2]
#> [1,] 0 1.000000
#> [2,] 1 1.289892
#> [3,] 2 1.968273
#> [4,] 3 2.703592
#> [5,] 4 2.899549
#> [6,] 5 3.880089
#> [7,] 6 4.621610
#> [8,] 7 4.673057
#> [9,] 8 5.203269
#> [10,] 9 5.899093
#> [11,] 10 6.587649
#>
#> , , 4
#>
#> [,1] [,2]
#> [1,] 0 1.000000
#> [2,] 1 1.031230
#> [3,] 2 1.256793
#> [4,] 3 1.557624
#> [5,] 4 2.194089
#> [6,] 5 2.673114
#> [7,] 6 3.105285
#> [8,] 7 3.811719
#> [9,] 8 4.760296
#> [10,] 9 4.940634
#> [11,] 10 5.157534
#>
#> , , 5
#>
#> [,1] [,2]
#> [1,] 0 1.000000
#> [2,] 1 1.680163
#> [3,] 2 2.179009
#> [4,] 3 2.820688
#> [5,] 4 3.480972
#> [6,] 5 3.576996
#> [7,] 6 4.342597
#> [8,] 7 5.112271
#> [9,] 8 6.102984
#> [10,] 9 7.073505
#> [11,] 10 7.462687
#>
#> , , 6
#>
#> [,1] [,2]
#> [1,] 0 1.000000
#> [2,] 1 1.461186
#> [3,] 2 1.776428
#> [4,] 3 1.951104
#> [5,] 4 2.482678
#> [6,] 5 2.976315
#> [7,] 6 3.755623
#> [8,] 7 3.959802
#> [9,] 8 4.673199
#> [10,] 9 4.738415
#> [11,] 10 5.092622
#>
#> , , 7
#>
#> [,1] [,2]
#> [1,] 0 1.000000
#> [2,] 1 1.825199
#> [3,] 2 2.099018
#> [4,] 3 2.669063
#> [5,] 4 3.004782
#> [6,] 5 3.601044
#> [7,] 6 3.792563
#> [8,] 7 4.740326
#> [9,] 8 5.282807
#> [10,] 9 5.827410
#> [11,] 10 6.106007
#>
#> , , 8
#>
#> [,1] [,2]
#> [1,] 0 1.000000
#> [2,] 1 1.446702
#> [3,] 2 1.818214
#> [4,] 3 1.846275
#> [5,] 4 2.312262
#> [6,] 5 2.702293
#> [7,] 6 2.722358
#> [8,] 7 3.099329
#> [9,] 8 3.659242
#> [10,] 9 4.516326
#> [11,] 10 4.901135
#>
#> , , 9
#>
#> [,1] [,2]
#> [1,] 0 1.000000
#> [2,] 1 1.527917
#> [3,] 2 2.128555
#> [4,] 3 2.389926
#> [5,] 4 2.679976
#> [6,] 5 3.160051
#> [7,] 6 4.080057
#> [8,] 7 4.480777
#> [9,] 8 4.693950
#> [10,] 9 5.365717
#> [11,] 10 5.424331
#>
#> , , 10
#>
#> [,1] [,2]
#> [1,] 0 1.000000
#> [2,] 1 1.997069
#> [3,] 2 2.146105
#> [4,] 3 2.664661
#> [5,] 4 3.510781
#> [6,] 5 4.229051
#> [7,] 6 4.470365
#> [8,] 7 5.017408
#> [9,] 8 5.852210
#> [10,] 9 5.880166
#> [11,] 10 6.349551
#>