[R-sig-dyn-mod] Individual Based Model type with deSolve

[Romain Pete]

Hi Folks,

I’m trying to implement an IBM type model in an existing biogeochemical model (NPZD-like). What I got is an embedded (NOT nested function) DEB oyster model (in my NPZD) on which I’d like to use a for loop to solve DEB equations for multiple seeding. Finally, I’d like to get the output of DEB with the output of the NPZD:

MyModel <- function(t, state, parameters) {
	with(as.list(c(state, parameters)), {

# all biogeochemical processes defined here

for( s in seeding){
	run DEB equations
	c(seeding1, seeding2, etc..) # 4 state variables DEB hence a matrix result with nrow=4 state variables, ncol= « s » seeding
}

# rate of change
dN = remineralisation + oyster excretion - uptake
dP = growth - Zoo_grazing - Oyster_grazing - death
dZ = growth - Oyster_grazing - death - egestion
dD = death + egestion - mineralization

dE = stock			###
dV = volume			###	These 4 state variables are multiplied by number of seedings
dEr = reproduction		###
dEgo = gamete		###

list(c(
	dN, dP, dZ, dD, dO, c(seeding)
	))
	})
}


Trouble here is that the ode function uses initialization vector of length = number of state variables and when using multiple seedings, I, somehow, add state variables, hence, I’m struggling…

Any thoughts would be much appreciated

Cheers

Romain



> Le 22 nov. 2017 à 10:39, Maciek Jacek Swat <maciej.swat at gmail.com> a écrit :
> 
> hi All,
> 
> I ran a simple ODE model (see PS for R and Matlab code) and get order of
> magnitude longer running times with lsoda (as implemented in deSolve) than
> with ode45 in Matlab. (lsoda takes 2.97 seconds, ode45 takes 0.12 seconds).
> Is this so or can it be optimized?
> 
> For example, the ODEs in R are coded using the notation dA <- ..., dB <-
> ... and not the vectorized one, is it the reason for the longer simulation
> times?
> 
> Any comments would be very appreciated,
> Best, Maciej
> 
> 
> PS:
> R code:
> rm(list=ls())
> 
> Nanda <- function(t, state, parameters) {
>  with(as.list(c(state, parameters)), {
> 
>    dB_CLL <- 1/body*(((r*B_CLL)*body) - ((d_B*B_CLL)*body) + ((b_B)*body)
> - ((d_BN*B_CLL*NK)*body) - ((d_BT*B_CLL*Tc)*body))
>    dNK <- 1/body*(-((d_NB*NK*B_CLL)*body) + ((b_N)*body) -
> ((d_N*NK)*body))
>    dTc <- 1/body*(-((d_TB*Tc*B_CLL)*body) - ((d_T*Tc)*body) + ((b_T)*body)
> + (k*a_TH*B_CLL^L/(s+B_CLL^L)*TH*Tc))
>    dTH <- 1/body*((a_TH*B_CLL^L/(s+B_CLL^L)*TH) - ((d_TH*TH)*body) +
> ((b_TH)*body))
> 
>    list(c(
>      dB_CLL,dNK,dTc,dTH
>    ))
>  })
> }
> 
> parameters <- c(
>  a_TH = 0.0001,
>  L = 2,
>  s = 1000,
>  d_B = 0.0076,
>  r = 0.0176,
>  b_B = 68,
>  B_CLL_0 = 100,
>  NK_0 = 1000,
>  Tc_0 = 100,
>  TH_0 = 100,
>  d_BN = 0.0001,
>  d_BT = 0.026,
>  d_NB = 0.0001,
>  b_N = 3.98,
>  d_N = 0.0159,
>  d_TB = 0.001,
>  d_T = 0.005,
>  b_T = 2.6,
>  k = 0.6,
>  d_TH = 0.01,
>  b_TH = 10.4,
>  body = 1
> )
> 
> state <- c(
>  B_CLL = 100,
>  NK = 1000,
>  Tc = 100,
>  TH = 100
> )
> 
> t1=proc.time()
> times <- seq(0, 300, by = 0.01)
> out <- ode(y = state, times = times, func = Nanda, method="lsoda", parms =
> parameters)
> cat((proc.time()-t1)[3],'\n')
> 
> 
> 
> 
> Matlab code:
> clear all; close all; clc;
> 
> yy0(1) = 100; % B_CLL
> yy0(2) = 1000; % NK
> yy0(3) = 100; % Tc
> yy0(4) = 100; % TH
> 
> tic
> [T,Y] = ode45(@(t,y)NandaFCT(t,y),[0 300],yy0);
> toc
> 
> figure
> semilogy(T,Y); grid on;
> 
> 
> function dydt = NandaFCT(t,y)
> 
> a_TH = 0.0001;
> L = 2;
> s = 1000;
> d_B = 0.0076;
> r = 0.0176;
> b_B = 68;
> B_CLL_0 = 100;
> NK_0 = 1000;
> Tc_0 = 100;
> TH_0 = 100;
> d_BN = 0.0001;
> d_BT = 0.026;
> d_NB = 0.0001;
> b_N = 3.98;
> d_N = 0.0159;
> d_TB = 0.001;
> d_T = 0.005;
> b_T = 2.6;
> k = 0.6;
> d_TH = 0.01;
> b_TH = 10.4;
> body = 1;
> 
> dydt = zeros(4,1);
> 
> dydt(1)= 1/body*(((r*y(1))*body) - ((d_B*y(1))*body) + ((b_B)*body) -
> ((d_BN*y(1)*y(2))*body) - ((d_BT*y(1)*y(3))*body)); % B_CLL
> dydt(2)= 1/body*(-((d_NB*y(2)*y(1))*body) + ((b_N)*body) -
> ((d_N*y(2))*body)); % NK
> dydt(3)= 1/body*(-((d_TB*y(3)*y(1))*body) - ((d_T*y(3))*body) +
> ((b_T)*body) + (k*a_TH*y(1)^L/(s+y(1)^L)*y(4)*y(3))); % Tc
> dydt(4)= 1/body*((a_TH*y(1)^L/(s+y(1)^L)*y(4)) - ((d_TH*y(4))*body) +
> ((b_TH)*body)); % TH
> 
> end
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models