[R-sig-dyn-mod] better implementation of variable species ode.2D model

[Jeremy Chacon]

Hello list,

I am running ode.2D models (from deSolve, using ReacTran to manage
diffusion) in which I wish to have a variable number of "species"
(differential equations) which a single function can handle.  I've managed
to come up with a method that uses lists and loops internally, see the code
below for an example with a simpler reaction than I will be using in
practice.

However, it's quite slow, especially when I scale up the number of species
(or the number of boxes, but that is expected).

I would greatly appreciate any tips on how to better implement this idea.
I know if it was a "mean-field" model, I could simply use matrix algebra
inside the function to allow a variable number of species, but I'm not sure
how that can be extended to a 2D situation.

Thanks for any tips!

Jeremy

Code:

# example of variable-species ode.2D model
rm(list = ls())
library(ReacTran)

variable_species_lv_2D = function(t, y, parms){


  # y holds all the data.  there are Nx * Ny boxes.
  n_colonies = length(y) / (Nx*Ny) - 1

  # the resource is always first in the data vector y
  resource = matrix(nrow = Nx, ncol = Ny,
              data = y[1:(Nx * Ny)])

  # grab each consumer's data
  consumers = list()
  for (k in 1:n_colonies){
    consumers[[k]] = matrix(nrow = Nx, ncol = Ny,
                    data = y[(1 + k * (Nx * Ny)): ((k + 1) * Nx * Ny)])
  }

  # get the diffusion for the resource
  resource_diff = tran.2D(resource, grid = grid2D, D.grid =
D_resource_grid)$dC

  # get the diffusion for each consumer
  consumers_diff = list()
  for (k in 1:n_colonies){
    consumers_diff[[k]] = tran.2D(consumers[[k]], grid = grid2D, D.grid =
D_consumer_grid)$dC
  }

  # get the reaction for each consumer, and simultaneously
  # subtract that from the resource
  consumers_rxn = list()
  resource_rxn = matrix(nrow = Nx, ncol = Ny,
                  data = 0)
  for (k in 1:n_colonies){
    consumers_rxn[[k]] = r * consumers[[k]] * resource
    resource_rxn = resource_rxn - consumers_rxn[[k]]
  }

  # sum the effects of diffusion + reaction,
  # and simultaneously concatenate all the data back
  # into a single vector
  dresource = resource_diff + resource_rxn
  dAll = dresource
  for (k in 1:n_colonies){
    curr_dconsumer = consumers_diff[[k]] + consumers_rxn[[k]]
    dAll = c(dAll, curr_dconsumer)
  }
  return(list(dAll))
}

# setup parameters
world_width = 100
Nx = Ny = 30
D_resource = 1
D_consumer = 1e-2

r = 0.01

# generate spatial vectors / matrices
Gridx = setup.grid.1D(L = world_width, N = Nx)
Gridy = setup.grid.1D(L = world_width, N = Ny)
grid2D = setup.grid.2D(Gridx, Gridy)
D_consumer_grid = setup.prop.2D(value = D_consumer, y.value = D_consumer,
grid = grid2D)
D_resource_grid = setup.prop.2D(value = D_resource, y.value = D_resource,
grid = grid2D)


# initialize data
n_consumers = 4
resource_start = matrix(nrow = Nx, ncol = Ny, data = 1000)
all_start = resource_start
for (k in 1:n_consumers){
  curr_consumer_start = matrix(nrow = Nx, ncol = Ny, data = 0)
  curr_consumer_start[sample(1:(Nx * Ny), 1)] = 1
  all_start = c(all_start, curr_consumer_start)
}

# run model
times = seq(from = 0, to = 10, by = 1)
out = ode.2D(y = all_start,
             times = times,
            func = variable_species_lv_2D,
            parms = NULL,
            nspec = n_consumers + 1,
            dimens = c(Nx, Ny),
            lrw = 5e8,
            hmax = 0.01)

# visualize output
# image(out, mfrow = c(1, n_consumers + 1))

-- 

*___________________________________________________________________________Jeremy
M. Chacon, Ph.D.*

*Post-Doctoral Associate, Harcombe Lab*
*University of Minnesota*
*Ecology, Evolution and Behavior*

	[[alternative HTML version deleted]]