Alternative to for loop for fast calculations when equations depend on each other

user_2829961 0
Viewed 116

I am using a for-loop to do step-by-step calculations where several equations depend on each other. Because of this dependence, I cannot find a solution where I do the calculations inside a dataframe. My main motivation is to speed up the calculations when the Time vector is very large in the reprex below.

Could you please suggest alternatives to the following for-loop based calculations, preferably inside a dataframe in R? The only thing I can think of is using for-loop in Rcpp.

Reproducible Example

last_time <- 10
STEP = 1

Time <- seq(from = 0, to = last_time, by = STEP)

## empty vectors 
eq1 <- vector(mode = "double", length = length(Time)) 
eq2 <- vector(mode = "double", length = length(Time)) 
eq <- vector(mode = "double", length = length(Time)) 
eq3 <- vector(mode = "double", length = length(Time)) 
eq4 <- vector(mode = "double", length = length(Time)) 

## adding the first values
eq1[1] <- 25
eq2[1] <- 25
eq[1] <- 25 
eq3[1] <- 100
eq4[1] <- 2

for (t in 2:length(Time)) { 
  ## eq1
  eq1[t] <- eq[t-1] + (2.5 *  STEP * (1 - (eq[t-1])/25))
  ## eq2
  eq2[t] <- (-2 * STEP) + ((-2^2) * (STEP^2)) - (2 * eq3[t-1]) - (eq[t-1] * STEP)
  ## min.
  eq[t] <- min(eq1[t], eq2[t] )
  ## eq3
  eq3[t] <- (eq[t] - eq[t-1])/(STEP)
  ## eq4
  eq4[t] <- eq4[t-1] + (eq[t-1] * STEP) + (0.5 * eq3[t-1] * (STEP)^2)
}

Output:

my_data <- data.frame(Time, eq1, eq2, eq, eq3, eq4) 
my_data
#>    Time        eq1        eq2         eq        eq3        eq4
#> 1     0   25.00000   25.00000   25.00000 -256.00000     2.0000
#> 2     1   25.00000 -231.00000 -231.00000   25.60000  -101.0000
#> 3     2 -205.40000  225.00000 -205.40000   23.04000  -319.2000
#> 4     3 -182.36000  199.40000 -182.36000   20.73600  -513.0800
#> 5     4 -161.62400  176.36000 -161.62400   18.66240  -685.0720
#> 6     5 -142.96160  155.62400 -142.96160   16.79616  -837.3648
#> 7     6 -126.16544  136.96160 -126.16544   15.11654  -971.9283
#> 8     7 -111.04890  120.16544 -111.04890   13.60489 -1090.5355
#> 9     8  -97.44401  105.04890  -97.44401   12.24440 -1194.7819
#> 10    9  -85.19961   91.44401  -85.19961   11.01996 -1286.1037
#> 11   10  -74.17965   79.19961  -74.17965    0.00000 -1365.7934

Created on 2021-02-28 by the reprex package (v1.0.0)

dataframe r
2 Answers

You could define a recursive function. A loop is faster than recursion though.

g <- function(m, STEP, time, x=2) {
  if (time == 0) m
  else {
    ## eq1
    m[x, 2] <- m[x - 1, 1] + 2.5*STEP*(1 - (m[x - 1, 1])/25)
    ## eq2
    m[x, 3] <- -2*STEP + -2^2*STEP^2 - 2*m[x - 1, 4] - m[x - 1, 1]*STEP
    ## min.
    m[x, 1] <- min(m[x, 2], m[x, 3])
    ## eq3
    m[x - 1, 4] <- (m[x, 1] - m[x - 1, 1])/STEP
    ## eq4
    m[x, 5] <- m[x - 1, 5] + m[x - 1, 1]*STEP + 0.5*m[x - 1, 4]*STEP^2
    g(m, STEP, time - 1, x + 1)
  }
}

Usage

last_time <- 10; STEP <- 1

First <- c(eq0=25, eq1=25, eq2=25, eq3=100, eq4=2)
m <- matrix(0, last_time + 1, length(First), dimnames=list(NULL, names(First)))
m[1, ] <- First
    
g(m, STEP, last_time)
#              eq0        eq1        eq2        eq3        eq4
#  [1,]   25.00000   25.00000   25.00000 -256.00000     2.0000
#  [2,] -231.00000   25.00000 -231.00000   25.60000  -101.0000
#  [3,] -205.40000 -205.40000  225.00000   23.04000  -319.2000
#  [4,] -182.36000 -182.36000  199.40000   20.73600  -513.0800
#  [5,] -161.62400 -161.62400  176.36000   18.66240  -685.0720
#  [6,] -142.96160 -142.96160  155.62400   16.79616  -837.3648
#  [7,] -126.16544 -126.16544  136.96160   15.11654  -971.9283
#  [8,] -111.04890 -111.04890  120.16544   13.60489 -1090.5355
#  [9,]  -97.44401  -97.44401  105.04890   12.24440 -1194.7819
# [10,]  -85.19961  -85.19961   91.44401   11.01996 -1286.1037
# [11,]  -74.17965  -74.17965   79.19961    0.00000 -1365.7934
user_65740380

as you asked how it works: The recursive filter function of stats::filter can be used with mapply as follows:

dataframe <-  
    mapply(stats::filter, 
           dataframe, 
           filter = vector, 
           method = "recursive")

where vector is e.g. c(25), which could be your first eq1[1] <- 25

The recursive filter works like a recursive loop but is a bit more elegant: Then the mapply recursive filter would do:

                   dataframe / vector
row or timepoint 1  20 
row or timepoint 2  30       + (20 * c(25))
row or timepoint 3  40       + ((20*25)+30) * c(25))

It calculates the value in the first row and uses it in the next, where it multiplies the next vector. Perhaps if you play around with stats filter and the recursive method you also get the same result. It is a row based calculation over time similar to Rcpp but more flexible.

user_139942620
Related