porting Matlab cells to R to collect output from FOR loop

Question

I have an ODE model in Matlab for which I'm interested in performing some parameter sweeps.

I am trying to port the following code from Matlab to R

for i = 1:numel(sweep1)
  initial_conditions(6)=sweep1(i);
    for j = 1:numel(sweep2) 
      parameters(3)=sweep2(j); 
[t,y] = ode23s(@(timespan, initial_conditions) MODEL(timespan, initial_conditions,          parameters), timespan, initial_conditions);
    results_cell{i,j}=[y(end,1),y(end,2)];

The 2 FOR statements above vary first 1 initial condition (i), then for each i vary a parameter (j) and run the solver. The output from the solver for each iteration of the loop is then collected in a cell 'results_cell'

This runs fine in Matlab but I need to port it to R. The loops are the same and the solver code is implemented using deSolve, however I am not sure how to collect the results from the solver at each iteration of the loop as R doesn't have cells like Matlab, and how to gather {i,j} from each loop along with the 2 ode outputs.

Ultimately I would like to plot a heat map of the ode solver output vs the values in each of the 2 parameter sweeps.

Thanks for any help.

Answer 1

I'm assuming sweep1 and sweep2 are both vectors of numbers. What you can do is use expand.grid to make a data frame of the combinations of that, and then loop over the frame once with apply:

# sweep 1, sweep 2
sweep1 <- c(1, 2, 4)
sweep2 <- c(3, 5, 7)

# expand out the combinations
combinations <- expand.grid(sweep1=sweep1, sweep2=sweep2)

# apply over the data frame
results <- apply(combinations, 1, function(row) {
  # set up the parameters from the row which has been passed in.
  initial_conditions[6] <- row["sweep1"]
  parameters[3]         <- row["sweep2"]

  # call ode23s 
  res <- ode23s(initial_conditons, parameters, function, whatever, ...)

  # there should be a nicer way than calling nrow twice here, but R doesn't
  # seem to have the nice 'end' keyword
  # also, we copy in the row, so that's in the output.
  c(row, one=res[nrow(res), 1], two=res[nrow(res), 2])
})

# because the apply has flipped rows to columns...
results <- as.data.frame(t(results))

results
#     sweep1 sweep2 one  two
#  1  1      3      ...  ...
#  2  2      3      ...  ...
# ...

The result of all this is a data frame of the input combinations and the output combinations. If you want more factors, add on a sweep3, but beware of the combinatorial complexity...

Answer 2

Here what I would do: I run the ode23 once to get the structure of the solution.

sweep1 =2
sweep2 =3
library(pracma)
f <- function(t, x,i=1,j=0)
  as.matrix(c(x[1] * ((i+j) - x[2]^2) -x[2], x[1]))
t0 <- 0
tf <- 20
x0 <- as.matrix(c(0, 0.25))
sol = ode23(f, t0, tf, x0,1,1)$y
res = tail(sol,1)

Then I use replicate to create the structure of the final output matrix. Using this trick avoid us to deal with pre-allocating arrays. replicate will do for us.

results_cell = replicate(sweep1,replicate(sweep2,res))

I just run my final simulation and assign each solution to results_cell

for (i in seq(sweep1))
  for (j in seq(sweep2))
    results_cell[,,j,i] = tail(ode23(f, t0, tf, x0,i,j)$y,1)