search

kingaa / pomp

R package for statistical inference using partially observed Markov processes
https://kingaa.github.io/pomp
GNU General Public License v3.0
110 stars 26 forks source link

Error: in ‘simulate’: Mistaken array variable declaration in pomp, maybe? #122

Closed luismieryteran closed 3 years ago

luismieryteran commented 3 years ago

Hi there,

I am new to the pomp package and to C, potentially deadly combinations. I'm trying to set up a disease transmission model with an explicit serial interval in it. To force it to be Markovian, I truncate the serial interval at 5 days (for now) and make the state variable an array of length 5. Each step of the Markov process outputs the new 5 days worth of incidence given the previous 5 days, etc.

The model is here:

# -----------------------------------------------------------
#       Parameters
# -----------------------------------------------------------
stateSize <- 5

nSteps <- round(360 / stateSize)
pop <- 10e6

R0 <- 2.5
car <- 0.05
phi_car <- 5

I0 <- 10

# -----------------------------------------------------------
#       Simulate
# -----------------------------------------------------------
rinit_C <- Csnippet("
  double *I = &I1;
  double *E_Cases = &E_Cases1;

  int d;
  int N = (int) stateSize;
  for (d = 0; d < N - 1; d++) {
    I[d] = 0.0;
    E_Cases[d] = 0.0;
  }
  I[N - 1] = I0;
  E_Cases[N - 1] = 0.0;

  cumsumI = 0.0;")

rproc_C <- Csnippet("
  int N = (int) stateSize;

  double serialIntervalDistributionNorm = 0.0;
  double serialIntervalDistribution[2 * N];

  double carDelayDistributionNorm = 0.0;
  double carDelayDistribution[2 * N];

  int d;
  for(d = 0; d < N; d++){
    serialIntervalDistribution[d] = 0.1;
    serialIntervalDistribution[d + N] = 0.0;
    serialIntervalDistributionNorm += serialIntervalDistribution[d];

    carDelayDistribution[d] = 0.1;
    carDelayDistribution[d + N] = 0.0;
    carDelayDistributionNorm += carDelayDistribution[d];
  }

  // Normalizing distributions
  for(d = 0; d < 2 * N; d++){
    serialIntervalDistribution[d] = serialIntervalDistribution[d] / serialIntervalDistributionNorm;

    carDelayDistribution[d] = carDelayDistribution[d] / carDelayDistributionNorm;
  }

  double *I = &I1;
  double *E_Cases = &E_Cases1;
  double I_temp[2 * N];
  double E_Cases_temp[2 * N];

  for (d = 0; d < N; d++) {
    I_temp[d] = I[d];
    E_Cases_temp[d] = E_Cases[d];
  }
  for (d = N; d < 2 * N; d++) {
    I_temp[d] = 0.0;
    E_Cases_temp[d] = 0.0;
  }

  int d1, d2;
  double newI;
  for (d1 = N;  d1 < 2 * N; d1++){

    for (d2 = 1; d2 <= N; d2++){
      newI = R0 * (1 - cumsumI / pop) * I_temp[d1 - d2] * serialIntervalDistribution[d2 - 1];
      if (newI < 0) newI = 0;

      I_temp[d1] += newI;
      cumsumI += newI;

      E_Cases_temp[d1] += car * I_temp[d1 - d2] * carDelayDistribution[d2 - 1];

    }
  }

  for (d1 = 0;  d1 < N; d1++){
    I[d1] = I_temp[N + d1];
    E_Cases[d1] = E_Cases_temp[d1];
  }")

rmeas_C <- Csnippet("
  int N = (int) stateSize;

  double *I = &I1;
  double *Cases = &Cases1;

  int d;
  for (d = 0; d < N; d++){
    Cases[d] = rnbinom(phi_car, phi_car / (phi_car + car * I[d]))
  }")

sim <- simulate(t0 = 0, 
                times = 1:nSteps,
                statenames = c(sprintf("I%d",seq_len(stateSize)), 
                               "cumsumI", 
                               sprintf("E_Cases%d",seq_len(stateSize))),
                obsnames = c(sprintf("Cases%d", seq_len(stateSize))),
                paramnames = c("pop", "R0", "car", "phi_car", "I0", "stateSize"),
                params = c(pop = pop, R0 = R0, car = car, phi_car = phi_car, I0 = I0,
                           stateSize = as.integer(stateSize)),
                rinit = rinit_C,
                rprocess = discrete_time(step.fun = rproc_C, delta.t = 1),
                rmeasure = rmeas_C,
                nsim = 100) %>% as.data.frame()

The error that results when trying to run the above is

Error: in ‘simulate’: error in building shared-object library from C snippets: in ‘Cbuilder’: compilation error: cannot compile shared-object library ‘/var/folders/pz/qwn1wrxn33s5rnqjqpqbmlth0000gn/T//Rtmpx277ts/94987/pomp_fd965b5ba09ff864df576a7c23ee4369.so’: status = 1
compiler messages:
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I/Library/Frameworks/R.framework/Versions/4.0/Resources/library/pomp/include -I/Users/luis/Dropbox/Documents/Learning/Particle_Filter  -I/usr/local/include   -fPIC  -Wall -g -O2  -c /var/folders/pz/qwn1wrxn33s5rnqjqpqbmlth0000gn/T//Rtmpx277ts/94987/pomp_fd965b5ba09ff864df576a7c23ee4369.c -o /var/folders/pz/qwn1wrxn33s5rnqjqpqbmlth0000gn/T//Rtmpx277ts/94987/pomp_fd965b5ba09ff864df576a7c23ee4369.o
/var/folders/pz/qwn1wrxn33s5rnqjqpqbmlth0000gn/T//Rtmpx277ts/94987/pomp_fd965b5ba09ff864df576a7c23ee4369.c:197:11: warning: initializing 'double *' with an expression of type 'const double *' discards
In addition: Warning message:
In system2(command = R.home("bin/R"), args = c("CMD", "SHLIB", "-c",  :
  running command 'PKG_CPPFLAGS="-I/Library/Frameworks/R.framework/Versions/4.0/Resources/library/pomp/include -I/Users/luis/Dropbox/Documents/Learning/Particle_Filter" '/Library/Frameworks/R.framework/Resources/bin/R' CMD SHLIB -c -o /var/folders/pz/qwn1wrxn33s5rnqjqpqbmlth0000gn/T//Rtmpx277ts/94987/pomp_fd965b5ba09ff864df576a7c23ee4369.so /var/folders/pz/qwn1wrxn33s5rnqjqpqbmlth0000gn/T//Rtmpx277ts/94987/pomp_fd965b5ba09ff864df576a7c23ee4369.c 2>&1' had status 1

and the output of source("https://kingaa.github.io/scripts/diagnostics.R") is

-----------------------------------
 contents of user Renviron file 
-----------------------------------

-----------------------------------
 sessionInfo 
-----------------------------------
R version 4.0.2 (2020-06-22)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.6

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] doRNG_1.8.2       rngtools_1.5      doParallel_1.0.15 iterators_1.0.12  foreach_1.5.0     subplex_1.6       pomp_3.1.1.4      tictoc_1.0       
 [9] forcats_0.5.0     dplyr_1.0.0       purrr_0.3.4       readr_1.3.1       tidyr_1.1.0       tibble_3.0.3      ggplot2_3.3.2     tidyverse_1.3.0  
[17] stringr_1.4.0    

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.5        lubridate_1.7.9   mvtnorm_1.1-1     lattice_0.20-41   assertthat_0.2.1  digest_0.6.25     utf8_1.1.4        R6_2.4.1         
 [9] cellranger_1.1.0  plyr_1.8.6        backports_1.1.8   reprex_0.3.0      coda_0.19-3       httr_1.4.2        pillar_1.4.6      rlang_0.4.7      
[17] readxl_1.3.1      rstudioapi_0.11   blob_1.2.1        labeling_0.3      munsell_0.5.0     broom_0.7.0       compiler_4.0.2    modelr_0.1.8     
[25] pkgconfig_2.0.3   tidyselect_1.1.0  codetools_0.2-16  viridisLite_0.3.0 fansi_0.4.1       crayon_1.3.4      dbplyr_1.4.4      withr_2.2.0      
[33] grid_4.0.2        jsonlite_1.7.0    gtable_0.3.0      lifecycle_0.2.0   DBI_1.1.0         magrittr_1.5      scales_1.1.1      cli_2.0.2        
[41] stringi_1.4.6     farver_2.0.3      reshape2_1.4.4    fs_1.4.2          xml2_1.3.2        ellipsis_0.3.1    generics_0.0.2    vctrs_0.3.2      
[49] deSolve_1.28      tools_4.0.2       glue_1.4.1        hms_0.5.3         colorspace_1.4-1  rvest_0.3.5       haven_2.3.1

Perhaps you guys can help me see what is amiss? Really appreciate the help.

kingaa commented 3 years ago

Thanks @luismieryteran, for such a carefully prepared question. There were just a few small errors in your C code. See the following:

library(tidyverse)
library(pomp)

rmeas_C <- Csnippet("
  int N = (int) stateSize;

  const double *I = &I1;
  double *Cases = &Cases1;

  int d;
  for (d = 0; d < N; d++){
    Cases[d] = rnbinom(phi_car, phi_car / (phi_car + car * I[d]));
  }")

simulate(
  t0 = 0,
  times = 1:nSteps,
  statenames = c(sprintf("I%d",seq_len(stateSize)),
                 "cumsumI",
                 sprintf("E_Cases%d",seq_len(stateSize))),
  obsnames = c(sprintf("Cases%d", seq_len(stateSize))),
  paramnames = c("pop", "R0", "car", "phi_car", "I0", "stateSize"),
  params = c(pop = pop, R0 = R0, car = car, phi_car = phi_car, I0 = I0,
             stateSize = stateSize),
  rinit = rinit_C,
  rprocess = discrete_time(step.fun = rproc_C, delta.t = 1),
  rmeasure = rmeas_C,
  nsim = 100
) %>%
  as.data.frame() -> sim

sim %>%
  select(time,rep=.id,starts_with("Cases")) %>%
  gather(var,cases,-time,-rep) %>%
  group_by(time,var) %>%
  summarize(
    p=c(0.05,0.5,0.95),
    q=quantile(cases,probs=p),
    lab=c("lo","med","hi")
  ) %>%
  ungroup() %>%
  select(-p) %>%
  spread(lab,q) %>%
  ggplot(aes(x=time,y=med,ymin=lo,ymax=hi))+
  geom_ribbon(fill=grey(0.8),color=NA)+
  geom_line()+
  guides(color=FALSE)+
  facet_grid(var~.)

This runs for me and should for you. A couple of pointers:

  1. Remember: he semicolon ; is needed to end a line in C.
  2. Read up on the const declaration in C and cf. your error messages.

Am closing this issue, but feel free to reopen if more discussion is warranted.