How to set limits for a parameter? #514
Replies: 2 comments 3 replies
-
|
Hi @Ricardo-DO, Thank you for reaching out. I modified your original post so the code in your comment is interpreted as R code so it looks a little nicer and the You suggested using a Depending on the version of library(rxode2)
#> Warning: package 'rxode2' was built under R version 4.2.3
#> rxode2 2.0.13 using 4 threads (see ?getRxThreads)
#> no cache: create with `rxCreateCache()`
f <- function() {
ini({
z = 0.008
V = 250
X0 = 20
Cs = 0.1
})
model({
d/dt(Xv) = z*X0*(((X0-Xv)/X0)^(1/3))*(Cs-(Xv/V))
})
}
ev <- eventTable()
ev$add.sampling(0:600) #in seconds
x <- rxSolve(f, ev)
# runs fine
print(x)
#> ── Solved rxode2 object ──
#> ── Parameters ($params): ──
#> z V X0 Cs
#> 0.008 250.000 20.000 0.100
#> ── Initial Conditions ($inits): ──
#> Xv
#> 0
#> ── First part of data (object): ──
#> # A tibble: 601 × 2
#> time Xv
#> <dbl> <dbl>
#> 1 0 0
#> 2 1 0.0160
#> 3 2 0.0320
#> 4 3 0.0479
#> 5 4 0.0639
#> 6 5 0.0798
#> # ℹ 595 more rows
# update to Cs = 1.6
f2 <- f |> ini(Cs=1.6)
#> ℹ change initial estimate of `Cs` to `1.6`
x <- rxSolve(f2, ev)
# runs fine, but some values are NaN
summary(x)
#> ── Summary of Solved rxode2 object ──
#> ── Model ($model): ──
#> rxode2({
#> param(z, V, X0, Cs)
#> d/dt(Xv) = z * X0 * (((X0 - Xv)/X0)^(1/3)) * (Cs - (Xv/V))
#> })
#> ── Parameters (x$params): ──
#> z V X0 Cs
#> 0.008 250.000 20.000 1.600
#> ── Initial Conditions (x$inits): ──
#> Xv
#> 0
#> ── Summary of data (object): ──
#> time Xv
#> Min. : 0 Min. : 0.00
#> 1st Qu.:150 1st Qu.: 7.05
#> Median :300 Median :12.91
#> Mean :300 Mean :11.95
#> 3rd Qu.:450 3rd Qu.:17.40
#> Max. :600 Max. :19.96
#> NA's :481
# replace the line in the mdel with ifelse()
f3 <- f2 |>
model({d/dt(Xv) <-ifelse(Xv<X0, z*X0*(((X0-Xv)/X0)^(1/3))*(Cs-(Xv/V)),X0)})
print(f3)
#> ── rxode2-based free-form 1-cmt ODE model ──────────────────────────────────────
#> ── Initalization: ──
#> Fixed Effects ($theta):
#> z V X0 Cs
#> 0.008 250.000 20.000 1.600
#>
#> States ($state or $stateDf):
#> Compartment Number Compartment Name
#> 1 1 Xv
#> ── Model (Normalized Syntax): ──
#> function() {
#> ini({
#> z <- 0.008
#> V <- 250
#> X0 <- 20
#> Cs <- 1.6
#> })
#> model({
#> d/dt(Xv) <- ifelse(Xv < X0, z * X0 * (((X0 - Xv)/X0)^(1/3)) *
#> (Cs - (Xv/V)), X0)
#> })
#> }
x <- rxSolve(f3, ev)
# No values are NaN by using ifelse
summary(x)
#> ── Summary of Solved rxode2 object ──
#> ── Model ($model): ──
#> rxode2({
#> param(z, V, X0, Cs)
#> d/dt(Xv) = ifelse(Xv < X0, z * X0 * (((X0 - Xv)/X0)^(1/3)) *
#> (Cs - (Xv/V)), X0)
#> })
#> ── Parameters (x$params): ──
#> z V X0 Cs
#> 0.008 250.000 20.000 1.600
#> ── Initial Conditions (x$inits): ──
#> Xv
#> 0
#> ── Summary of data (object): ──
#> time Xv
#> Min. : 0 Min. : 0.0
#> 1st Qu.:150 1st Qu.: 603.1
#> Median :300 Median :3603.1
#> Mean :300 Mean :3846.5
#> 3rd Qu.:450 3rd Qu.:6603.1
#> Max. :600 Max. :9603.1Created on 2023-04-29 with reprex v2.0.2 Note that the way you use library(rxode2)
#> Warning: package 'rxode2' was built under R version 4.2.3
#> rxode2 2.0.13 using 4 threads (see ?getRxThreads)
#> no cache: create with `rxCreateCache()`
#############################################
# Define model
ode <- "
d/dt(Xv) = z*X0*(((X0-Xv)/X0)^(1/3))*(Cs-(Xv/V));
"
# Compile model
mod1 <- RxODE(model = ode, modName = "mod1")
# Define system parameters
params <-
c(z = 0.008,
V = 250,
X0 = 20 ,
Cs = 0.1)
inits <- 0
# Initialize event table
ev <- eventTable()
# Specify sampling
ev$add.sampling(0:600) #in seconds
# Simulate
x <- mod1$run(params, ev, inits)
#> Warning: Assumed order of inputs: Xv
#####################################Created on 2023-04-29 with reprex v2.0.2 |
Beta Was this translation helpful? Give feedback.
-
|
Hi @mattfidler Many thanks for taking the time to reply to my question, really appreciate it! Thanks for tidying up the code, it's neater now. I get the same summary: ── Parameters (x$params): ── As we want that when Xv >= X0, the equation changes to X0, the max value should be X0, which in this case is 20. The plot looks like this
If I ran this similar code in Berkeley Madonna METHOD RK4 init xdv = 0 x0 = 20; d/dt (xdv) = zf*x0 * (((x0-xdv)/x0)^(1/3)) * (Cs - (xdv/Vv)) I get this plot
Do you think that there is something that it wasn't captured? Any advice would be greatly appreciated. Many thanks, |
Beta Was this translation helpful? Give feedback.
-
|
Sorry for the 2 replies. I kept trying things, and I made this change to the code (which I acknwledge is not the best solution, but the result is similar to what I'm looking for) If in the ifelse statement we keep the same equation but instead of X0-Xv we directly write 0 as the else condition, it kind of works library(rxode2) f <- function() { ev <- eventTable() ev$add.sampling(0:3600) #in seconds x <- rxSolve(f, ev)
── Summary of Solved rxode2 object ── We can see that now, the max should be 20 and it is 22.08, but it's closer to 20 in this way. Thanks, |
Beta Was this translation helpful? Give feedback.
-
|
This is probably what you are looking for: library(rxode2)
#> rxode2 2.0.13.9000 using 8 threads (see ?getRxThreads)
#> no cache: create with `rxCreateCache()`
f <- function() {
ini({
z = 0.008
V = 250
X0 = 20
Cs = 1.6
})
model({
d/dt(Xv) <- ifelse(Xv <= X0, z * X0 * (((X0 - Xv)/X0)^(1/3)) * (Cs - (Xv/V)), 0)
})
}
ev <- eventTable()
ev$add.sampling(0:3600) #in seconds
x <- rxSolve(f, ev)
library(ggplot2)
plot(x) + xlim(0, 600)
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
#> Warning: Removed 3000 rows containing missing values (`geom_line()`).Created on 2023-05-01 with reprex v2.0.2 |
Beta Was this translation helpful? Give feedback.
-
|
That's great @mattfidler! That is what I was lloking for! Best, |
Beta Was this translation helpful? Give feedback.
-
Hello everyone,
I was trying to use RxODE (by now, I am planning to move to rxode2 but I having problems with the installation) to replicate a code I was using in Berkeley Madonna.
It is a dissolution model that I am planning to implement in a compartmental model.
The code is the following
The code runs ok in those conditions, but if I change Cs to 1.6, Xv becomes larger than X0 after 100 s and I only get NaN as output.
Is there a way to specify an ifelse function like
Xv<-ifelse(Xv<X0, Xv,X0) ?
The problem is that, in those conditions X0-Xv will be <0
I tried to create another parameter, let's call it Xs = X0-Xv, but it says that Xs is needed to solve the system.
I tried creating the ifelse statement like
if (Xv<X0){
Xv=Xv} else{
Xv=X0}
But that doesn't seem to be the answer.
In Berkeley Madonna thee is a function called LIMIT which does that, is there something similar in RxODE?
If you have any thoughts it will be much appreciated.
Many thanks,
Ricardo
Beta Was this translation helpful? Give feedback.
All reactions