[R-sig-dyn-mod] Accessing internal time-varying and fixed parameters in root function and event function in deSolve

[Karline Soetaert]

Hi Martin, if I understand well, your root function does not directly depend on the state variables, but on another variable.

There are two things that you can do
1.	Recalculate the other variable (as a function of the state variables that are at your disposal) in the root function.

2.	Make the other variable a GLOBAL variable. This can be done by 
   a.	Define the variable first outside of the model function
   b.	Within the function, assign the new value with <<-  (which accesses the variable at a scope one higher)
   c.	Within the root function, access the global varialbe
For instance:
gE <- NULL # create the varlable globally

Seaweed_model <- function (,,,){
...
gE <<- ... # change the global variable

....
}

Rootfunc <- function(){...

...gE  # will access the global variable
}

This option 2 will work, but is not very neat

Option 1 is the most clean one, and I would implement it as such, i.e. recalculate the variable within the root function

All the best!


Karline Soetaert

-----Original Message-----
From: R-sig-dynamic-models <r-sig-dynamic-models-bounces using r-project.org> On Behalf Of Dr.Martin Johnson
Sent: woensdag 26 januari 2022 23:28
To: r-sig-dynamic-models using r-project.org
Subject: [R-sig-dyn-mod] Accessing internal time-varying and fixed parameters in root function and event function in deSolve

[You don't often get email from mjohnson using bmrs.ie. Learn why this is important at http://aka.ms/LearnAboutSenderIdentification.]

Hi all,

I'm using R+desolve to model seaweed growth and wish to implement a function in which seaweed is harvested given a threshold value of an internal variable of the model (i.e. not a state variable so not accessed via the 'y' of (t, y, parms)). I would also like the threshold value to be controlled by a variable that I can pass to the model in parms.  Some (simplified) example code

Seaweed_model <- function(t, y, parms,...) {

  with(as.list(c(y, parms)), {
    #(lots of internal variables hidden here for brevity)

   I_top<-PAR(t)*exp(-K_d*(z-h_MA))                                                          # calculate incident irradiance at the top of the farm
    I_av <-(I_top/(K_d*z+N_f*a_cs))*(1-exp(-(K_d*z+N_f*a_cs)))          # calclulate the average irradiance throughout the height of the farm
    g_E <- I_av/((I_s)+I_av))                                                                          # light limitation scaling function


#(state variables not relevant to this problem omitted)

    dN_f        <- mu_g_EQT*N_s-(d_m*N_f)                                    # change in fixed nitrogen

    list(c(dNH4, dNO3,dN_s,dN_f,dD),
         g_E       = g_E,
        )
  })
}

Then I would like to trigger an event function when g_E reaches a threshold value

rootfunc  <- function(t,y,parms,...){
  return(g_E - threshold_value)  #how do I access these values from the ode solver call?
}

And the event

eventfunc <- function(t, y, parms,...){
  c(y[1:2],y[3]*harvest_fraction,y[4]*harvest_fraction,y[5]) #reduces biomass state variables fractionally }

The harvest parameters are added in to the parms passed to the ode call

parms_harvest<-c(
  harvest_threshold=0.2,
  harvest_fraction=0.75
)

lim_harvest <- ode(times = times, func = Seaweed_model,  y = y0, parms = c(parms_porphyra,parms_farm,parms_harvest), events=list(func=eventfunc, root=TRUE),rootfun=rootfunc)

Am I trying to achieve the impossible? All of the examples for desolve events do a root evaluation on a state variable of the model (e.g. y['N_s'] is accessible in the example above) rather than internal variables of the model. g_E as written in the above rootfunc clearly will not work but I have no idea how to access it if it's even possible.  Substituting  a state variable of the model for g_E (not useful scientifically but good for testing purposes) I get error: object 'threshold_value' not found so it doesn't seem like the root function is able to access parms even.

I'm sure I've missed something fundamental here. Thanks so much in advance for you help!

Martin Johnson




        [[alternative HTML version deleted]]

_______________________________________________
R-sig-dynamic-models mailing list
R-sig-dynamic-models using r-project.org
https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-dynamic-models&data=04%7C01%7C%7Cda75e0c551c24720986c08d9e164b2d5%7C9a1651bf58af435b86a83e9334b4b732%7C0%7C0%7C637788645349264000%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=m0GFOj6l5QSCPCUqTuC2do%2BmJ6kzXsEgB%2FfCvITHCoM%3D&reserved=0