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

[Jeremy Chacon]

Thank you Thomas for your reply (not to mention your fantastic packages!).

I tried changing my function so that the data all get tranformed between a
1d vector and a (preallocated) 3d array (new code below), but when I
compared the two head-to-head, found no difference between my original
loop-list attempt and the newer loop-array attempt. So perhaps I am simply
hitting the point where, as you suggest, I move the function itself outside
of R.  That said, I still get the sense I am doing something fundamentally
slow with my code, but can't quite figure out what that is.

Thank you again,

Jeremy

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

# list implementation (old way)
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))
}


# matrix implementation (new way)

diffusion_matrix = array(dim = c(Nx, Ny, n_consumers+1), data = 0)
reaction_matrix = array(dim = c(Nx, Ny, n_consumers+1), data = 0)
array_lv_2D = function(t, y, parms){

  all_data_array = array(y, dim = c(Nx, Ny, n_consumers + 1))

  # get the diffusion for the resource
  diffusion_matrix[,,1] = tran.2D(all_data_array[,,1], grid = grid2D,
D.grid = D_resource_grid)$dC

  # get the diffusion for each consumer
  for (k in 2:(n_consumers+1)){
    diffusion_matrix[,,k] = tran.2D(all_data_array[,,k], grid = grid2D,
D.grid = D_consumer_grid)$dC
  }

  # get the reaction for each consumer, and simultaneously
  # subtract that from the resource
  for (k in 2:(n_consumers+1)){
    reaction_matrix[,,k] = r * all_data_array[,,k] * all_data_array[,,1]
    reaction_matrix[,,1] = reaction_matrix[,,1] - reaction_matrix[,,k]
  }
  return(list(as.vector(diffusion_matrix + reaction_matrix)))
}



world_width = 100
Nx = Ny = 40
D_resource = 1
D_consumer = 1e-2

r = 0.01

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)


n_consumers = 8
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)
}

times = seq(from = 0, to = 10, by = 1)

print(system.time(
  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)
))

print(system.time(
  out <- ode.2D(y = all_start,
                times = times,
                func = array_lv_2D,
                parms = NULL,
                nspec = n_consumers + 1,
                dimens = c(Nx, Ny),
                lrw = 5e8,
                hmax = 0.01)
))

# note both times are the same


On Fri, Feb 24, 2017 at 10:43 AM, Thomas Petzoldt <
thomas.petzoldt at tu-dresden.de> wrote:

> Am 24.02.2017 um 16:31 schrieb 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.
>>
>
> 2D works nicely if you flatten your 2D matrix into a 1D vector (with
> "as.vector") and back to a matrix with "matrix" when needed. Reshaping
> matrices is very fast and my experience shows that also deSolve with
> multidimensional models can be fast, if vectorization is consequently used.
> If this is not possible, you may consider to use compiled code (C or
> Fortran).
>
> Regards, Thomas
>
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
>



-- 

*___________________________________________________________________________Jeremy
M. Chacon, Ph.D.*

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

	[[alternative HTML version deleted]]