-xde_solve(mosy1,Tmax=50,dt=10)$outputs$orbits$deout -> out2xde_solve(mosy1,Tmax=50,dt=10)$outputs$orbits$deout -> out2
+#> Called from: get_bionomics_s(tm, de_vars, pars, 1)
+#> debug: pars <- reset_state_i(1, tm, de_vars, pars)
+#> debug: f <- get_f(pars, s)
+#> debug: q <- get_q(pars, s)
+#> debug: g <- get_g(pars, s)
+#> debug: sigma <- get_sigma(pars, s)
+#> debug: Omega <- compute_Omega(pars, s)
+#> debug: Upsilon <- compute_Upsilon(pars, s)
+#> debug: for (i in 2:length(tm)) {
+#> pars <- reset_state_i(i, tm, de_vars, pars)
+#> f <- rbind(f, get_f(pars, s))
+#> q <- rbind(q, get_q(pars, s))
+#> g <- rbind(g, get_g(pars, s))
+#> sigma <- rbind(sigma, get_sigma(pars, s))
+#> }
+#> debug: pars <- reset_state_i(i, tm, de_vars, pars)
+#> debug: f <- rbind(f, get_f(pars, s))
+#> debug: q <- rbind(q, get_q(pars, s))
+#> debug: g <- rbind(g, get_g(pars, s))
+#> debug: sigma <- rbind(sigma, get_sigma(pars, s))
+#> debug: pars <- reset_state_i(i, tm, de_vars, pars)
+#> debug: f <- rbind(f, get_f(pars, s))
+#> debug: q <- rbind(q, get_q(pars, s))
+#> debug: g <- rbind(g, get_g(pars, s))
+#> debug: sigma <- rbind(sigma, get_sigma(pars, s))
+#> debug: pars <- reset_state_i(i, tm, de_vars, pars)
+#> debug: f <- rbind(f, get_f(pars, s))
+#> debug: q <- rbind(q, get_q(pars, s))
+#> debug: g <- rbind(g, get_g(pars, s))
+#> debug: sigma <- rbind(sigma, get_sigma(pars, s))
+#> debug: pars <- reset_state_i(i, tm, de_vars, pars)
+#> debug: f <- rbind(f, get_f(pars, s))
+#> debug: q <- rbind(q, get_q(pars, s))
+#> debug: g <- rbind(g, get_g(pars, s))
+#> debug: sigma <- rbind(sigma, get_sigma(pars, s))
+#> debug: pars <- reset_state_i(i, tm, de_vars, pars)
+#> debug: f <- rbind(f, get_f(pars, s))
+#> debug: q <- rbind(q, get_q(pars, s))
+#> debug: g <- rbind(g, get_g(pars, s))
+#> debug: sigma <- rbind(sigma, get_sigma(pars, s))
+#> debug: return(list(f = f, q = q, g = g, sigma = sigma))
diff --git a/docs/articles/ex_534.html b/docs/articles/ex_534.html
index 78da4ee2b..ecc5a97bd 100644
--- a/docs/articles/ex_534.html
+++ b/docs/articles/ex_534.html
@@ -184,6 +184,7 @@ 5-3-4 Example
library(data.table) library(ggplot2) library(deSolve) +#> ℹ Loading ramp.xdsIntroduction
@@ -225,31 +226,31 @@Encoding Structural Information
For the aquatic habitats, the habitat membership vector is an ordered
list of the index of patches where the habitats are found:
-
+
membership = c(1,1,1,2,2)
membership
#> [1] 1 1 1 2 2
The number of habitats, \(n_q\) or
nHabitats, is the length of the membership matrix:
-
+
nHabitats = length(membership)
nHabitats
#> [1] 5
For the human population strata, the residence vector is an ordered
list of the index of patches where people live:
-
+
residence = c(2,2,3,3)
residence
#> [1] 2 2 3 3
The number of strata, \(n_h\) or
nStrata, is the length of the residence matrix:
-
+
nStrata = length(residence)
nStrata
#> [1] 4
The number of patches, \(n_p\) or
nPatches, is just a number:
-
+
nPatches = 3
@@ -268,7 +269,7 @@ Aquatic Habitat Membership Matrix
The matrix is created by a function
create_habitat_matrix:
-
+
habitat_matrix = create_habitat_matrix(nPatches, membership)
habitat_matrix
#> [,1] [,2] [,3] [,4] [,5]
@@ -277,7 +278,7 @@ Aquatic Habitat Membership Matrix#> [3,] 0 0 0 0 0
The habitat matrix is used to sum quantities in habitats up to
patches. For example, the number of habitats per patch is \({\cal N} \cdot 1\), where \(1\) is a vector of 1’s:
-
@@ -297,7 +298,7 @@ The Residence Matrix
+
residence_matrix = create_residence_matrix(nPatches, residence)
residence_matrix
#> [,1] [,2] [,3] [,4]
@@ -305,7 +306,7 @@ The Residence Matrix#> [2,] 1 1 0 0
#> [3,] 0 0 1 1
The number of strata per patch is:
-
@@ -321,14 +322,14 @@ Egg Laying\(Q\).
For now, we generate arbitary weights for each one of the
habitats:
-
+
searchWtsQ = c(7,2,1,8,2)
searchWtsQ
#> [1] 7 2 1 8 2
And we can compute availability as \(\cal N
\cdot w\). In ramp.xds, the function that computes
habitat availability is called compute_Q:
-
+
Q <- compute_Q(habitat_matrix, searchWtsQ)
Q
#> [1] 10 10 0
@@ -355,7 +356,7 @@ Egg Laying
+
calU = compute_calU(searchWtsQ, habitat_matrix, Q)
calU
#> [,1] [,2] [,3]
@@ -380,10 +381,10 @@ Blood Feeding
+
searchWtsH = c(1,1,1,1)
This size of each population stratum differs:
-
+
HPop <- c(10,90, 100, 900)
The time spent matrix, \(\Theta\) is
a \(n_p \times n_h\) matrix. Each
@@ -399,7 +400,7 @@
Blood Feeding
+
TaR <- t(matrix(
c(c(0.01,0.01,0.001,0.001),
c(.95,.92,.04,.02),
@@ -407,7 +408,7 @@ Blood Feeding))
The availability of hosts is of time at risk and the population
density weighed by their search weights:
-
+
W <- compute_W(searchWtsH, HPop, TaR)
W
#> [1] 2.0 114.3 934.2
@@ -418,7 +419,7 @@ Transmission\(\beta\),
describes the expected proportion of each infective bite that would be
received be each population stratum.
-
+
compute_beta(HPop, W, searchWtsH, TaR)
#> [,1] [,2] [,3]
#> [1,] 5e-03 0.0083114611 4.281738e-05
@@ -439,7 +440,7 @@ Aquatic Mosquito DynamicsLpar
.
-
+
psi <- rep(1/8, nHabitats)
phi <- rep(1/8, nHabitats)
theta <- c(1/10, 1/20, 1/40, 1/100, 1/10)
@@ -456,7 +457,7 @@ Adult Mosquito Dynamics
The parameter values are:
-
+
g <- 1/12
sigma <- 1/12/2
mu <- 0
@@ -466,12 +467,12 @@ Adult Mosquito Dynamicseip <- 12
eggsPerBatch <- 30
We create a named list:
-
+
MYZo = list(g=g, sigma=sigma, mu=mu, f=f, q=q, nu=nu, eip=eip, eggsPerBatch=eggsPerBatch)
Mosquito dispersal among the patches is described by a matrix, \(\cal K.\) Each column in \(\cal K\) describes the proportion of
emigrating mosquitoes that go to every other patch. The diagonal
elements are all \(0\):
-
+
calK <- t(matrix(
c(c(0, .6, .3),
c(.4, 0, .7),
@@ -480,21 +481,21 @@ Adult Mosquito Dynamics It is computed by
compute_Omega_xde:
-
+
Omega <- compute_Omega_xde(g, sigma, mu, calK)
Survival and dispersal through the EIP in this model, denoted \(\Upsilon,\) is computed using matrix
exponentiation:
-
+
Upsilon <- expm::expm(-Omega*eip)
The function xde_make_MYZpar_RM can be used to construct
the adult mosquito model object, called MYZpar.
-
+
MYZpar = create_MYZpar_RM(nPatches, MYZo)
class(MYZpar) <- "RM"
The parameters are assigned to a list called baseline so
that it can be stored and used to compute the values of bionomic
parameters that have been modified by control.
-
+
names(MYZpar)
#> [1] "nPatches" "eip_par" "eip" "f_t" "q_t"
#> [6] "g_t" "sigma_t" "mu" "nu" "eggsPerBatch"
@@ -509,12 +510,12 @@ Human Infection DynamicsThe 5-3-4 model uses the basic SIS (Susceptible-Infected-Susceptible)
model for the human component (see more here). It can be configured using
ramp.xds::xde_make_Xpar_SIS.
-
+
r <- 1/200
b <- 0.55
c <- c(0.1, .02, .1, .02)
The model is configured and assigned the name Xpar:
-
+
Xpar <- create_Xpar_SIS(nStrata, list(), b, r, c)
@@ -522,12 +523,12 @@ Initial Conditions
+
For the adult mosquito model:
-
+
MYZinits = list(
M = rep(100, nPatches),
P = rep(10, nPatches),
@@ -535,7 +536,7 @@ Initial Conditions= rep(0, nPatches)
)
For the human model:
-
@@ -609,15 +610,15 @@ The xds Templatemembership and residence, as explained
above.
-
+
params = make_xds_template("xde", "full", "dde", nPatches, membership, residence)
After being set up:
-
+
c(params$nHabitats, params$nPatches, params$nStrata)
#> [1] 5 3 4
This was created by make_xds_template and stored as
params$habitat_matrix
-
+
params$habitat_matrix
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 1 1 1 0 0
@@ -625,7 +626,7 @@ The xds Template#> [3,] 0 0 0 0 0
If we want to retrieve the membership matrix, we can call
view_habitat_matrix
-
+
view_habitat_matrix(params)
#> $habitat_index
#> [1] 1 2 3 4 5
@@ -636,7 +637,7 @@ The xds TemplateQ
and calU get created with
default values: Iall patches are assumed to have the same biting
weights.
-
+
params$vars$Q
#> [[1]]
#> [1] 3 2 0
@@ -645,7 +646,7 @@ The xds Templatemake_xds_template only sets up the first species.
Similarly, we can view the egg distribution matrix for the first
species:
-
+
params$calU[[1]]
#> [,1] [,2] [,3]
#> [1,] 0.3333333 0.0 0
@@ -666,23 +667,23 @@ The xds Template
Building the Object
-
+
params$Lpar = list()
params$Lpar[[1]] = Lpar
params$MYZpar = list()
params$MYZpar[[1]] = MYZpar
params$Xpar = list()
params$Xpar[[1]] = Xpar
-
+
params <- make_Linits(params, 1, Linits)
params <- make_MYZinits(params, 1, MYZinits)
params <- make_Xinits(params, HPop, 1, Xinits)
-
+
params <- setup_Hpar_static(params, 1)
After the parameters for 5-3-4 have been specified, we
can generate the indices for the model and attach them to the parameter
list.
-
+
params = make_indices(params)
@@ -693,15 +694,15 @@ Setting Parameter ValuesmakeQ and make_calU do this.
-
+
params <- change_habitat_weights(params, searchWtsQ)
params <- make_Q(params)
params <- make_calU(params)
We can check to see
-
+
params$vars$Q[[1]]
#> [1] 10 10 0
-
+
params$calU[[1]]
#> [,1] [,2] [,3]
#> [1,] 0.7 0.0 0
@@ -709,17 +710,17 @@ Setting Parameter Values#> [3,] 0.1 0.0 0
#> [4,] 0.0 0.8 0
#> [5,] 0.0 0.2 0
-
+
params <- change_TimeSpent(TaR, params)
params <- make_TaR(params)
-
+
params$TimeSpent
#> [[1]]
#> [,1] [,2] [,3] [,4]
#> [1,] 0.01 0.01 0.001 0.001
#> [2,] 0.95 0.92 0.040 0.020
#> [3,] 0.04 0.02 0.959 0.929
-
+
params$TaR
#> [[1]]
#> [[1]][[1]]
@@ -727,27 +728,27 @@ Setting Parameter Values#> [1,] 0.01 0.01 0.001 0.001
#> [2,] 0.95 0.92 0.040 0.020
#> [3,] 0.04 0.02 0.959 0.929
-
+
params <- change_blood_weights(params, searchWtsH)
y0 <- get_inits(params)
params <- make_WB(0, params, y0)
-
+
params$vars$W[[1]]
#> [1] 2.0 114.3 934.2
-
+
params <- change_calK(calK, params)
params$MYZpar[[1]]$calK
#> [,1] [,2] [,3]
#> [1,] 0.0 0.6 0.3
#> [2,] 0.4 0.0 0.7
#> [3,] 0.6 0.4 0.0
-
+
get_Omega(params, 1)
#> [,1] [,2] [,3]
#> [1,] 0.12500000 -0.02500000 -0.01250000
#> [2,] -0.01666667 0.12500000 -0.02916667
#> [3,] -0.02500000 -0.01666667 0.12500000
-
+
get_Upsilon(params, 1)
#> [,1] [,2] [,3]
#> [1,] 0.23643097 0.07241259 0.04640076
@@ -780,7 +781,7 @@ Setting Parameter ValuesSolving
Now we can get the initial conditions of the model.
-
+
params <- make_indices(params)
@@ -795,14 +796,14 @@ Numerical Solution
+
We want to pass an unnamed vector to the solver so:
-
-y0 = get_inits(params, flatten=TRUE)
+y0 = get_inits(params, flatten=TRUE)
+
out <- deSolve::dede(y = y0, times = seq(0, 735, by =15),
func = xde_derivatives, parms = params)
out1 <- out
@@ -812,7 +813,7 @@ Plot Output
+
colnames(out)[params$ix$L[[1]]$L_ix+1] <- paste0('L_', 1:params$nHabitats)
colnames(out)[params$ix$MYZ[[1]]$M_ix+1] <- paste0('M_', 1:params$nPatches)
colnames(out)[params$ix$MYZ[[1]]$P_ix+1] <- paste0('P_', 1:params$nPatches)
@@ -837,7 +838,7 @@ Plot OutputUsing xde_setup
We create lists with all our parameters values:
-
+
MYZo = list(
solve_as = "ode",
g = 1/12, sigma = 1/12/2,
@@ -846,20 +847,20 @@ Using xde_setup eip = 12,
M = 100, P = 10, Y = 1, Z = 0
)
-
+
-
+
-
+
xds_setup(dlay='dde',
MYZname="RM", Xname="SIS", Lname="basicL",
nPatches = 3, HPop=c(10, 90, 100, 900),
@@ -868,17 +869,17 @@ Using xde_setup residence=c(2,2,3,3), searchB=searchWtsH,
TimeSpent=TaR, searchQ = c(7,2,1,8,2)) -> mod534
We solve and take the differences to check:
-
+
mod534 <- xde_solve(mod534, Tmax=735, dt=15)
mod534$outputs$orbits$deout -> out2
-
+
xds_plot_M(mod534, llty = c(1:3))
Interestingly, the differences are small:
-
+
-
+
approx_equal(tail(out2, 1), tail(out1,1), tol = 1e-5)
#> logical(0)
diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml
index 1032dba9d..bc6f9171e 100644
--- a/docs/pkgdown.yml
+++ b/docs/pkgdown.yml
@@ -26,7 +26,7 @@ articles:
Understanding_ramp.xds: Understanding_ramp.xds.html
vc_lemenach: vc_lemenach.html
Working: Working.html
-last_built: 2024-08-26T23:27Z
+last_built: 2024-08-29T19:50Z
urls:
reference: https://dd-harp.github.io/ramp.xds/reference
article: https://dd-harp.github.io/ramp.xds/articles
diff --git a/docs/reference/index.html b/docs/reference/index.html
index 6927c43cc..8b5452b66 100644
--- a/docs/reference/index.html
+++ b/docs/reference/index.html
@@ -232,6 +232,10 @@ Solve
xde_stable_orbit()
Solve for the steady state or stable orbit of a system of equations
+
+
+ Cohort orbits for a model \(\cal X\)
diff --git a/vignettes/GettingStarted.Rmd b/vignettes/GettingStarted.Rmd
index 498e4dbec..cea25ec6a 100644
--- a/vignettes/GettingStarted.Rmd
+++ b/vignettes/GettingStarted.Rmd
@@ -53,13 +53,14 @@ Note that `model` was passed to `xde_solve,` and the return value replaces `mode
names(model$outputs)
```
-The outputs are:
+Standard outputs are:
+ `time` -- are the time points when the values of the dependent variables were evaluated
+ `last_y` -- is the final state
+ `orbits` -- are the solutions with the variables parsed by name
+ `terms` -- the dynamical terms describing transmission, including the EIR.
+ `deout` -- the raw (unparsed) outputs
+
+ `bionomics` -- For models with exogenous forcing, `xde_solve` stores the values of mosquito bionomic parameters
When `xde_solve` is called, the **`xds` model object** it returns has replaced any old values. To do anything with the outputs, the user will probably have to write a wrapper function that solves and extracts.
+membership = c(1,1,1,2,2) membership #> [1] 1 1 1 2 2The number of habitats, \(n_q\) or
-nHabitats, is the length of the membership matrix:+nHabitats = length(membership) nHabitats #> [1] 5For the human population strata, the residence vector is an ordered list of the index of patches where people live:
-+residence = c(2,2,3,3) residence #> [1] 2 2 3 3The number of strata, \(n_h\) or
-nStrata, is the length of the residence matrix:+nStrata = length(residence) nStrata #> [1] 4The number of patches, \(n_p\) or
-nPatches, is just a number:+nPatches = 3@@ -268,7 +269,7 @@Aquatic Habitat Membership Matrix
The matrix is created by a function
-create_habitat_matrix:+habitat_matrix = create_habitat_matrix(nPatches, membership) habitat_matrix #> [,1] [,2] [,3] [,4] [,5] @@ -277,7 +278,7 @@Aquatic Habitat Membership Matrix#> [3,] 0 0 0 0 0
The habitat matrix is used to sum quantities in habitats up to patches. For example, the number of habitats per patch is \({\cal N} \cdot 1\), where \(1\) is a vector of 1’s:
- @@ -297,7 +298,7 @@The Residence Matrix
+residence_matrix = create_residence_matrix(nPatches, residence) residence_matrix #> [,1] [,2] [,3] [,4] @@ -305,7 +306,7 @@The Residence Matrix#> [2,] 1 1 0 0 #> [3,] 0 0 1 1
The number of strata per patch is:
- @@ -321,14 +322,14 @@Egg Laying\(Q\).
For now, we generate arbitary weights for each one of the habitats:
-+searchWtsQ = c(7,2,1,8,2) searchWtsQ #> [1] 7 2 1 8 2And we can compute availability as \(\cal N \cdot w\). In
-ramp.xds, the function that computes habitat availability is calledcompute_Q:+@@ -355,7 +356,7 @@Q <- compute_Q(habitat_matrix, searchWtsQ) Q #> [1] 10 10 0Egg Laying
+calU = compute_calU(searchWtsQ, habitat_matrix, Q) calU #> [,1] [,2] [,3] @@ -380,10 +381,10 @@Blood Feeding
+searchWtsH = c(1,1,1,1)This size of each population stratum differs:
-+HPop <- c(10,90, 100, 900)The time spent matrix, \(\Theta\) is a \(n_p \times n_h\) matrix. Each @@ -399,7 +400,7 @@
Blood Feeding
+TaR <- t(matrix( c(c(0.01,0.01,0.001,0.001), c(.95,.92,.04,.02), @@ -407,7 +408,7 @@Blood Feeding))
The availability of hosts is of time at risk and the population density weighed by their search weights:
-+@@ -418,7 +419,7 @@W <- compute_W(searchWtsH, HPop, TaR) W #> [1] 2.0 114.3 934.2Transmission\(\beta\), describes the expected proportion of each infective bite that would be received be each population stratum. -
+compute_beta(HPop, W, searchWtsH, TaR) #> [,1] [,2] [,3] #> [1,] 5e-03 0.0083114611 4.281738e-05 @@ -439,7 +440,7 @@. -Aquatic Mosquito DynamicsLpar
+psi <- rep(1/8, nHabitats) phi <- rep(1/8, nHabitats) theta <- c(1/10, 1/20, 1/40, 1/100, 1/10) @@ -456,7 +457,7 @@Adult Mosquito Dynamics
The parameter values are:
-+g <- 1/12 sigma <- 1/12/2 mu <- 0 @@ -466,12 +467,12 @@Adult Mosquito Dynamicseip <- 12 eggsPerBatch <- 30
We create a named list:
-+MYZo = list(g=g, sigma=sigma, mu=mu, f=f, q=q, nu=nu, eip=eip, eggsPerBatch=eggsPerBatch)Mosquito dispersal among the patches is described by a matrix, \(\cal K.\) Each column in \(\cal K\) describes the proportion of emigrating mosquitoes that go to every other patch. The diagonal elements are all \(0\):
-+calK <- t(matrix( c(c(0, .6, .3), c(.4, 0, .7), @@ -480,21 +481,21 @@Adult Mosquito Dynamics It is computed by
compute_Omega_xde: -+Omega <- compute_Omega_xde(g, sigma, mu, calK)Survival and dispersal through the EIP in this model, denoted \(\Upsilon,\) is computed using matrix exponentiation:
-+Upsilon <- expm::expm(-Omega*eip)The function
-xde_make_MYZpar_RMcan be used to construct the adult mosquito model object, calledMYZpar.+MYZpar = create_MYZpar_RM(nPatches, MYZo) class(MYZpar) <- "RM"The parameters are assigned to a list called
-baselineso that it can be stored and used to compute the values of bionomic parameters that have been modified by control.+names(MYZpar) #> [1] "nPatches" "eip_par" "eip" "f_t" "q_t" #> [6] "g_t" "sigma_t" "mu" "nu" "eggsPerBatch" @@ -509,12 +510,12 @@Human Infection DynamicsThe 5-3-4 model uses the basic SIS (Susceptible-Infected-Susceptible) model for the human component (see more here). It can be configured using
ramp.xds::xde_make_Xpar_SIS.-+r <- 1/200 b <- 0.55 c <- c(0.1, .02, .1, .02)The model is configured and assigned the name
-Xpar:+Xpar <- create_Xpar_SIS(nStrata, list(), b, r, c)@@ -522,12 +523,12 @@Initial Conditions
+For the adult mosquito model:
-+MYZinits = list( M = rep(100, nPatches), P = rep(10, nPatches), @@ -535,7 +536,7 @@Initial Conditions= rep(0, nPatches) )
For the human model:
-@@ -609,15 +610,15 @@The
xdsTemplatemembership andresidence, as explained above. -+params = make_xds_template("xde", "full", "dde", nPatches, membership, residence)After being set up:
-+c(params$nHabitats, params$nPatches, params$nStrata) #> [1] 5 3 4This was created by
-make_xds_templateand stored asparams$habitat_matrix+params$habitat_matrix #> [,1] [,2] [,3] [,4] [,5] #> [1,] 1 1 1 0 0 @@ -625,7 +626,7 @@The
xdsTemplate#> [3,] 0 0 0 0 0If we want to retrieve the membership matrix, we can call
-view_habitat_matrix+view_habitat_matrix(params) #> $habitat_index #> [1] 1 2 3 4 5 @@ -636,7 +637,7 @@andThe
xdsTemplateQcalUget created with default values: Iall patches are assumed to have the same biting weights. -+@@ -645,7 +646,7 @@params$vars$Q #> [[1]] #> [1] 3 2 0The
xdsTemplatemake_xds_template only sets up the first species. Similarly, we can view the egg distribution matrix for the first species: -+params$calU[[1]] #> [,1] [,2] [,3] #> [1,] 0.3333333 0.0 0 @@ -666,23 +667,23 @@The
xdsTemplateBuilding the Object
-+-params$Lpar = list() params$Lpar[[1]] = Lpar params$MYZpar = list() params$MYZpar[[1]] = MYZpar params$Xpar = list() params$Xpar[[1]] = Xpar+-params <- make_Linits(params, 1, Linits) params <- make_MYZinits(params, 1, MYZinits) params <- make_Xinits(params, HPop, 1, Xinits)+params <- setup_Hpar_static(params, 1)After the parameters for
-5-3-4have been specified, we can generate the indices for the model and attach them to the parameter list.+params = make_indices(params)@@ -693,15 +694,15 @@Setting Parameter ValuesmakeQ and
make_calUdo this. -+params <- change_habitat_weights(params, searchWtsQ) params <- make_Q(params) params <- make_calU(params)We can check to see
-+-params$vars$Q[[1]] #> [1] 10 10 0+-params$calU[[1]] #> [,1] [,2] [,3] #> [1,] 0.7 0.0 0 @@ -709,17 +710,17 @@Setting Parameter Values#> [3,] 0.1 0.0 0 #> [4,] 0.0 0.8 0 #> [5,] 0.0 0.2 0
+-params <- change_TimeSpent(TaR, params) params <- make_TaR(params)+-params$TimeSpent #> [[1]] #> [,1] [,2] [,3] [,4] #> [1,] 0.01 0.01 0.001 0.001 #> [2,] 0.95 0.92 0.040 0.020 #> [3,] 0.04 0.02 0.959 0.929+-params$TaR #> [[1]] #> [[1]][[1]] @@ -727,27 +728,27 @@Setting Parameter Values#> [1,] 0.01 0.01 0.001 0.001 #> [2,] 0.95 0.92 0.040 0.020 #> [3,] 0.04 0.02 0.959 0.929
+-params <- change_blood_weights(params, searchWtsH) y0 <- get_inits(params) params <- make_WB(0, params, y0)+-params$vars$W[[1]] #> [1] 2.0 114.3 934.2+-params <- change_calK(calK, params) params$MYZpar[[1]]$calK #> [,1] [,2] [,3] #> [1,] 0.0 0.6 0.3 #> [2,] 0.4 0.0 0.7 #> [3,] 0.6 0.4 0.0+-get_Omega(params, 1) #> [,1] [,2] [,3] #> [1,] 0.12500000 -0.02500000 -0.01250000 #> [2,] -0.01666667 0.12500000 -0.02916667 #> [3,] -0.02500000 -0.01666667 0.12500000+get_Upsilon(params, 1) #> [,1] [,2] [,3] #> [1,] 0.23643097 0.07241259 0.04640076 @@ -780,7 +781,7 @@Setting Parameter ValuesSolving
Now we can get the initial conditions of the model.
-+params <- make_indices(params)@@ -795,14 +796,14 @@Numerical Solution
+We want to pass an unnamed vector to the solver so:
--y0 = get_inits(params, flatten=TRUE)++y0 = get_inits(params, flatten=TRUE)@@ -812,7 +813,7 @@out <- deSolve::dede(y = y0, times = seq(0, 735, by =15), func = xde_derivatives, parms = params) out1 <- outPlot Output
+colnames(out)[params$ix$L[[1]]$L_ix+1] <- paste0('L_', 1:params$nHabitats) colnames(out)[params$ix$MYZ[[1]]$M_ix+1] <- paste0('M_', 1:params$nPatches) colnames(out)[params$ix$MYZ[[1]]$P_ix+1] <- paste0('P_', 1:params$nPatches) @@ -837,7 +838,7 @@Plot OutputUsing
xde_setupWe create lists with all our parameters values:
-+-MYZo = list( solve_as = "ode", g = 1/12, sigma = 1/12/2, @@ -846,20 +847,20 @@Using
xde_setupeip = 12, M = 100, P = 10, Y = 1, Z = 0 )+ -+ -+xds_setup(dlay='dde', MYZname="RM", Xname="SIS", Lname="basicL", nPatches = 3, HPop=c(10, 90, 100, 900), @@ -868,17 +869,17 @@Using
xde_setupresidence=c(2,2,3,3), searchB=searchWtsH, TimeSpent=TaR, searchQ = c(7,2,1,8,2)) -> mod534We solve and take the differences to check:
-+-mod534 <- xde_solve(mod534, Tmax=735, dt=15) mod534$outputs$orbits$deout -> out2+xds_plot_M(mod534, llty = c(1:3))
Interestingly, the differences are small:
-+ -diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml index 1032dba9d..bc6f9171e 100644 --- a/docs/pkgdown.yml +++ b/docs/pkgdown.yml @@ -26,7 +26,7 @@ articles: Understanding_ramp.xds: Understanding_ramp.xds.html vc_lemenach: vc_lemenach.html Working: Working.html -last_built: 2024-08-26T23:27Z +last_built: 2024-08-29T19:50Z urls: reference: https://dd-harp.github.io/ramp.xds/reference article: https://dd-harp.github.io/ramp.xds/articles diff --git a/docs/reference/index.html b/docs/reference/index.html index 6927c43cc..8b5452b66 100644 --- a/docs/reference/index.html +++ b/docs/reference/index.html @@ -232,6 +232,10 @@+approx_equal(tail(out2, 1), tail(out1,1), tol = 1e-5) #> logical(0)Solve
xde_stable_orbit()+ Solve for the steady state or stable orbit of a system of equations
+ + +Cohort orbits for a model \(\cal X\)
diff --git a/vignettes/GettingStarted.Rmd b/vignettes/GettingStarted.Rmd index 498e4dbec..cea25ec6a 100644 --- a/vignettes/GettingStarted.Rmd +++ b/vignettes/GettingStarted.Rmd @@ -53,13 +53,14 @@ Note that `model` was passed to `xde_solve,` and the return value replaces `mode names(model$outputs) ``` -The outputs are: +Standard outputs are: + `time` -- are the time points when the values of the dependent variables were evaluated + `last_y` -- is the final state + `orbits` -- are the solutions with the variables parsed by name + `terms` -- the dynamical terms describing transmission, including the EIR. + `deout` -- the raw (unparsed) outputs + + `bionomics` -- For models with exogenous forcing, `xde_solve` stores the values of mosquito bionomic parameters When `xde_solve` is called, the **`xds` model object** it returns has replaced any old values. To do anything with the outputs, the user will probably have to write a wrapper function that solves and extracts.