canmod / macpan2

Rebuilding https://github.com/mac-theobio/McMasterPandemic/
https://canmod.github.io/macpan2/
GNU General Public License v3.0
2 stars 0 forks source link

updates to `convolution` engine function #197

Open jfree-man opened 3 months ago

jfree-man commented 3 months ago
  1. macpan2::convolution computes different values than manual convolution computation,
> library(macpan2)
> library(dplyr)
> # kernel
> k=c(0.5, 0.25, 0.25)
> X=1
> Y=0
> expr = list(
+     X ~ X + X/2
+   , Y ~ convolution(x=X,k=k)
+ )
> sim_data = simple_sims(
+   iteration_exprs = expr
+   , time_steps = 10
+   , mats = list(X = X, Y = Y, k = k)
+ ) |> arrange(matrix)
> # manual computation
> # i = 1 to 3 (up to 3 lags)
> X_sim = (sim_data 
+          %>% filter(matrix=="X")
+          %>% mutate(
+             lag_1 = lag(value,default = 0)
+           , lag_2 = lag(value,2,default = 0)
+           , lag_3 = lag(value,3,default = 0))
+ )
> manual_convolution = X_sim$value + k[1]*X_sim$lag_1 + k[2]*X_sim$lag_2 + k[3]*X_sim$lag_3
> macpan2_convolution = sim_data %>% filter(matrix=="Y") %>% pull(value)
> ## view differences
> cbind(manual_convolution, macpan2_convolution)
      manual_convolution macpan2_convolution
 [1,]            1.50000            1.000000
 [2,]            3.00000            1.750000
 [3,]            4.87500            2.625000
 [4,]            7.68750            3.937500
 [5,]           11.53125            5.906250
 [6,]           17.29688            8.859375
 [7,]           25.94531           13.289062
 [8,]           38.91797           19.933594
 [9,]           58.37695           29.900391
[10,]           87.56543           44.850586
  1. Segmentation fault occurs when the kernel k is left out,
    # expr = list(
    #     X ~ X + X/2
    #   , Y ~ convolution(X)
    # )
    # 
    # simple_sims(
    #   iteration_exprs = expr
    #   , time_steps = 10
    #   , mats = list(X = 1, Y = 0)
    # ) 
stevencarlislewalker commented 3 months ago

Thanks @jfree-man. I have one observation. If I modify one of the lines above to this:

manual_convolution = k[1] * X_sim$value + k[2] * X_sim$lag_1 + k[3] * X_sim$lag_2

I get a better agreement, although there is still an issue with how the first few outputs are treated:

      manual_convolution macpan2_convolution
 [1,]           0.750000            1.000000
 [2,]           1.500000            1.750000
 [3,]           2.625000            2.625000
 [4,]           3.937500            3.937500
 [5,]           5.906250            5.906250
 [6,]           8.859375            8.859375
 [7,]          13.289062           13.289062
 [8,]          19.933594           19.933594
 [9,]          29.900391           29.900391
[10,]          44.850586           44.850586

It looks like both our interpretations are consistent with the equation given in the ms:

image

I get my interpretation if i starts at 0 and yours if i starts at 1.

jfree-man commented 3 months ago

@stevencarlislewalker I was assuming i started at 0, and $\phi(0)=1$ so the first term is just $X(t)$.

If i starts at 1, I think it should be this,

manual_convolution = k[1]*X_sim$lag_1 + k[2]*X_sim$lag_2 + k[3]*X_sim$lag_3

Then we get a better agreement than my initial assumption, but it still needs some fixing.

      manual_convolution macpan2_convolution
 [1,]           0.000000            1.000000
 [2,]           0.750000            1.750000
 [3,]           1.500000            2.625000
 [4,]           2.625000            3.937500
 [5,]           3.937500            5.906250
 [6,]           5.906250            8.859375
 [7,]           8.859375           13.289062
 [8,]          13.289062           19.933594
 [9,]          19.933594           29.900391
[10,]          29.900391           44.850586