Skip to content

[R-sig-dyn-mod] forcing variables (dummy ?) problem

5 messages · Soetaert, Karline, Thomas Petzoldt, Fred Jean

#
Dear list
I'm testing the use of the simecol package to get involved in 
mechanistic modelling in biology with R. (I'm an octave user for 
differential equations and R user for statistics, but I'd like to shift 
to get everything in R)

Since one week (I'm a stubborn dummy) I'm trying to understand how 
things work with forcing (external) data, but I'm getting nut because I 
don't understand how to get it work when having forcing data /at each 
timestep/.
As a sort of minimal example, here is an attempt to get such forcing 
inspired from the lv3 model example :

##============================================
require(simecol)

lvfred <- odeModel(
   main = function(time, init, parms, inputs) {

     # calcul par interpolation de la valeur de s.in ? chaque valeur de
     # time
     s.in <- approxTime1(inputs, time, rule = 2)["s.in"]
     #cat(length(s.in)) # = 1 ? chaque time

     with(as.list(c(init, parms)),{
         ds <- s.in  - b*s*p + g*k
         dp <- c*s*p - d*k*p
         dk <- e*p*k - f*k
         list(c(ds, dp, dk), s.in)
       })
   },

   parms = c(b = 0.1, c = 0.1, d = 0.1, e = 0.1, f = 0.1, g = 0.05),
   times  = seq(from = 0, to = 365, by = 1),
   #
   # On peut rajouter des sous equations au mod?le dans un fichier avec
   # un nom du type "equations.file.R"
   #
   # equations = equations.file,
   #
   #
   inputs = as.matrix(
     data.frame(
       #time = c(0,   99, 100,  150, 200, 365),
       #s.in = c(0.1, 0.2, 0.2, 0.2, 0.2, 0.1)
       time = times,
       s.in = c(rep(0.1,99), rep(0.2,101), rep(0.1, 166))
       )
     ),
   init = c(s = 1, p = 1, k = 1), # substrate, producer, consumer
   solver = "lsoda"
   )

this code works perfectly but I don't understand how and where in the 
code I should  include e.g. a forcing temperature function depending on 
the vector /times/ as in the short code her after :

#
Tm <- 12 	# Temp?rature moyenne
Ampl <- 6 	# Amplitude de variation
Per <- 365	# P?riode de la sinuso?de
Faz <- 365/12	# Phase de la sinuso?de

tempC <-  Tm - Ampl * sin((2*pi*times / Per) - Faz)
#

And once I get this /tempC/ signal, how to include it ? Should I use the 
/equations slot/ ?

Usin tempC, I tried to get a /t.in/ forcing variable build with 
approxTime1() as for /s.in/ but didn't succeed so far... I tried to find 
examples with such forcings in simecolModels package, but none of them 
seemed to have such tricks.

Could you please tell me where to find solutions to those problems ?

Best whishes

Fred J.
#
Fred,

I you type in R:
library(deSolve)
?forcings

A help page will be opened showing how to use forcing functions in the deSolve package. deSolve is the differential equation solver package used in simecol.


Karline

-----Original Message-----
From: r-sig-dynamic-models-bounces at r-project.org [mailto:r-sig-dynamic-models-bounces at r-project.org] On Behalf Of Fred Jean
Sent: 17 October 2011 22:23
To: R-sig-dynamic-models at r-project.org
Subject: [R-sig-dyn-mod] forcing variables (dummy ?) problem

Dear list
I'm testing the use of the simecol package to get involved in mechanistic modelling in biology with R. (I'm an octave user for differential equations and R user for statistics, but I'd like to shift to get everything in R)

Since one week (I'm a stubborn dummy) I'm trying to understand how things work with forcing (external) data, but I'm getting nut because I don't understand how to get it work when having forcing data /at each timestep/.
As a sort of minimal example, here is an attempt to get such forcing inspired from the lv3 model example :

##============================================
require(simecol)

lvfred <- odeModel(
   main = function(time, init, parms, inputs) {

     # calcul par interpolation de la valeur de s.in ? chaque valeur de
     # time
     s.in <- approxTime1(inputs, time, rule = 2)["s.in"]
     #cat(length(s.in)) # = 1 ? chaque time

     with(as.list(c(init, parms)),{
         ds <- s.in  - b*s*p + g*k
         dp <- c*s*p - d*k*p
         dk <- e*p*k - f*k
         list(c(ds, dp, dk), s.in)
       })
   },

   parms = c(b = 0.1, c = 0.1, d = 0.1, e = 0.1, f = 0.1, g = 0.05),
   times  = seq(from = 0, to = 365, by = 1),
   #
   # On peut rajouter des sous equations au mod?le dans un fichier avec
   # un nom du type "equations.file.R"
   #
   # equations = equations.file,
   #
   #
   inputs = as.matrix(
     data.frame(
       #time = c(0,   99, 100,  150, 200, 365),
       #s.in = c(0.1, 0.2, 0.2, 0.2, 0.2, 0.1)
       time = times,
       s.in = c(rep(0.1,99), rep(0.2,101), rep(0.1, 166))
       )
     ),
   init = c(s = 1, p = 1, k = 1), # substrate, producer, consumer
   solver = "lsoda"
   )

this code works perfectly but I don't understand how and where in the code I should  include e.g. a forcing temperature function depending on the vector /times/ as in the short code her after :

#
Tm <- 12 	# Temp?rature moyenne
Ampl <- 6 	# Amplitude de variation
Per <- 365	# P?riode de la sinuso?de
Faz <- 365/12	# Phase de la sinuso?de

tempC <-  Tm - Ampl * sin((2*pi*times / Per) - Faz)
#

And once I get this /tempC/ signal, how to include it ? Should I use the 
/equations slot/ ?

Usin tempC, I tried to get a /t.in/ forcing variable build with 
approxTime1() as for /s.in/ but didn't succeed so far... I tried to find 
examples with such forcings in simecolModels package, but none of them 
seemed to have such tricks.

Could you please tell me where to find solutions to those problems ?

Best whishes

Fred J.
2 days later
#
Dear Fred,

during object creation, the values of the slots are not accessible by 
each other for principle reasons and this is the case where a special 
initialization mecahism "initfunc" comes into play. The "initfunc" gets 
the raw incomplete version of the object (called object in creation) and 
one has all methods to access and manipulate the slots.

See package vignette or http://www.jstatsoft.org/v22/i09 for an example.

For the "lvfred" example, different possibilities exist (see below), but 
for your other problem you don't need "initfunc", just add the 
parameters of the sine function to the parms vector:

parms = c( ......, Tm = 12, Ampl = 6, Per = 365, Faz <- 365/12, ....)

and call the signal in the "main" model function:

with(as.list(c(init, parms)),{
   tempC <-  Tm - Ampl * sin((2*pi*times / Per) - Faz)
   .......
})

For more complicated cases with real data, use something like the 
following, where lvfred2 with approxfun is more flexible and faster.

Hope it helps

Thomas P.


require(simecol)

lvfred <- odeModel(
   main = function(time, init, parms, inputs) {
     s.in <- approxTime1(inputs, time, rule = 2)["s.in"]
     with(as.list(c(init, parms)),{
         ds <- s.in  - b*s*p + g*k
         dp <- c*s*p - d*k*p
         dk <- e*p*k - f*k
         list(c(ds, dp, dk), s.in)
       })
   },

   parms = c(b = 0.1, c = 0.1, d = 0.1, e = 0.1, f = 0.1, g = 0.05),

   times  = seq(from = 0, to = 365, by = 1),
   init = c(s = 1, p = 1, k = 1), # substrate, producer, consumer
   solver = "lsoda",
   initfunc = function(obj) {
     tmp <- as.matrix(data.frame(
       time = times(obj),
       s.in = c(rep(0.1,99), rep(0.2,101), rep(0.1, 166))
     ))
     inputs(obj) <- tmp
     return(obj)
   }
)

##==============================================================================


lvfred2 <- odeModel(
   main = function(time, init, parms, inputs) {
     s.in <- inputs$signal(time)
     with(as.list(c(init, parms)),{
         ds <- s.in  - b*s*p + g*k
         dp <- c*s*p - d*k*p
         dk <- e*p*k - f*k
         list(c(ds, dp, dk), s.in = s.in)
       })
   },

   parms = c(b = 0.1, c = 0.1, d = 0.1, e = 0.1, f = 0.1, g = 0.05),

   times  = seq(from = 0, to = 365, by = 1),
   init = c(s = 1, p = 1, k = 1), # substrate, producer, consumer
   solver = "lsoda",
   initfunc = function(obj) {
     tmp <- approxfun(times(obj),
                     c(rep(0.1,99), rep(0.2,101), rep(0.1, 166)))
     inputs(obj) <- list(signal = tmp)
     obj
   }
)

system.time(plot(sim(lvfred)))


system.time(plot(sim(lvfred2)))
#
And the last example can be even simpler (and faster) when approxfun 
gets only the required values:

lvfred3 <- lvfred2

initfunc(lvfred3) <- function(obj) {
   tmp <- approxfun(c(0, 98, 199, max(times(obj))),
                    c(0.1, 0.1, 0.2, 0.1), method = "constant", f = 1)
   inputs(obj) <- list(signal = tmp)
   obj
}


system.time(plot(sim(lvfred3)))


Thomas P.
#
Dear Thomas and Karline

Many thanks for your avices and for giving attention to my newbie's 
problems.

I didn't progress a lot those last days (had many courses).
At this point, before those last two emails, I followed the advice of 
Karline and here is a piece of code. As you will see, I commented the 
code from what I understood. I will now try to take into account advice 
from Thomas.

Many thanks again to both of you

Fred

##============================================
require(simecol)

rm(list=ls())
# external parameters

# time for constructing functions and for length of simulations
t1  = seq(from = 0, to = 800, by = 1)

# temperature forcing function parameters
Tm <- 12       # mean temperature
TAmpl <- 6     # temperature variation amplitude
TPer <- 365    # temperature sine period
TFaz <- 365/6  # temperature sine phase

#substrate forcing
Sm <- 0.5      # mean substrate input
SAmpl <- 0.5   # substrate input amplitude
SPer <- 365/2  # substrate input sine period
SFaz <- 365/12 # substrate input sine phase

#### The model
lvfred <- odeModel(
   main = function(time, init, parms, inputs) {

     # interpolate forcings - output of those is one interpolated value
     # at each integration step
     s.in <- approxTime1(inputs, time, rule = 2)["s.in"]
     t.in <- approxTime1(inputs, time, rule = 2)["t.in"]

     # model equations
     with(as.list(c(init, parms)),{
         ds <- s.in  - t.in*b*s*p + g*k
         dp <- c*s*p - d*k*p
         dk <- e*p*k - f*k
         list(c(ds, dp, dk), s.in, t.in)
       })
   },

   # values of parameters
   parms = c(b = 0.1, c = 0.1, d = 0.1, e = 0.1, f = 0.1, g = 0.05),

   # duration of simulations and integration timestep
   times  = c(from = 0, to = 800, by = 0.1),

   # One may add sub-equations in a file
   # "equations.file.R" -- see daphnia example in simecolModel
   # equations = equations.file,

   inputs = as.matrix(
     data.frame(
       #time = c(0,   99, 100,  150, 200, 365),
       #s.in = c(0.1, 0.2, 0.2, 0.2, 0.2, 0.1)
       time = t1,
       s.in = Sm - SAmpl * sin((2*pi*t1 / SPer) - SFaz),
       t.in = Tm - TAmpl * sin((2*pi*t1 / TPer) - TFaz)
       )
     ),

   init = c(s = 1, p = 1, k = 1), # substrate, producer, consumer

   solver = "rk4"
   )
##============================================