search

mrc-ide / odin

ᚩ A DSL for describing and solving differential equations in R
https://mrc-ide.github.io/odin
Other
104 stars 13 forks source link

modify SIR model by 'odin.dust::odin_dust' version2.0 #304

Closed whzhuscu closed 1 year ago

whzhuscu commented 1 year ago

Thanks. I have the most recent versions:

But the code also can not run.

Sorry, I made a mistake in pressing the issue (modify SIR model by 'odin.dust::odin_dust' version) close. https://github.com/mrc-ide/odin/issues/303 https://github.com/mrc-ide/odin/issues?q=is%3Aissue+is%3Aclosed

whzhuscu commented 1 year ago

When I run the code:

gen <- odin::odin({
  ylag <- delay(y, tau)
  initial(y) <- 0.5
  deriv(y) <- 0.2 * ylag * 1 / (1 + ylag^10) - 0.1 * y
  tau <- user(10)
  output(ylag) <- ylag
})

The error is:

load：pkgbuild
Generating model in c
ℹ Re-compiling odin58a9525b (debug build)
── R CMD INSTALL ────────────────────────────────────────────────────────────────────────
─  installing *source* package 'odin58a9525b' ...
   ** using staged installation
   ** libs
   make: "C:\Users\朱\AppData\Local\Temp\RtmpGksAEM\file2c244ae362f7": No such file or directory
   make: *** No rule to make target '"C:\Users\朱\AppData\Local\Temp\RtmpGksAEM\file2c244ae362f7"'.  Stop.
   ERROR: compilation failed for package 'odin58a9525b'
─  removing 'C:/Users/朱/AppData/Local/Temp/RtmpGksAEM/devtools_install_2c244583273/odin58a9525b'
Error in `(function (command = NULL, args = character(), error_on_status = TRUE, …`:
! System command 'Rcmd.exe' failed
---
Exit status: 1
stdout & stderr: <printed>
---
Type .Last.error to see the more details.
richfitz commented 1 year ago

The username 朱 is most likely confusing the compiler. See the help for ?tempdir for how to set this:

 By default, ‘tmpdir’ will be the directory given by ‘tempdir()’.
 This will be a subdirectory of the per-session temporary directory
 found by the following rule when the R session is started.  The
 environment variables ‘TMPDIR’, ‘TMP’ and ‘TEMP’ are checked in
 turn and the first found which points to a writable directory is
 used: if none succeeds ‘/tmp’ is used.  The path must not contain
 spaces.  Note that setting any of these environment variables in
 the R session has no effect on ‘tempdir()’: the per-session
 temporary directory is created before the interpreter is started.
whzhuscu commented 1 year ago

Thanks, I changed a computer to run the code:

sir <- odin.dust::odin_dust({
  update(S) <- S - n_SI
  update(I) <- I + n_SI - n_IR
  update(R) <- R + n_IR

  n_IR <- rbinom(I, 1 - exp(-beta * I / (S + I + R)))
  n_SI <- rbinom(S, 1 - exp(-gamma))

  initial(S) <- S_ini
  initial(I) <- I_ini
  initial(R) <- 0

  S_ini <- user(1000)
  I_ini <- user(10)
  beta <- user(0.2)
  gamma <- user(0.1)
})

But there is error:

ℹ 20 functions decorated with [[cpp11::register]]
✔ generated file 'cpp11.R'
✔ generated file 'cpp11.cpp'
Re-compiling odin26987586

─  installing *source* package 'odin26987586' ... (440ms)

   ** using staged installation

   ** libs

   "D:/Program Files/R/rtools/rtools40/mingw64/bin/"g++ -std=gnu++11  -I"D:/PROGRA~1/R/R-42~1.1/include" -DNDEBUG  -I'D:/Program Files/R/R-4.2.1/library/cpp11/include'     -IC:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include -DHAVE_INLINE -fopenmp    -O2 -Wall  -mfpmath=sse -msse2 -mstackrealign -c cpp11.cpp -o cpp11.o

   "D:/Program Files/R/rtools/rtools40/mingw64/bin/"g++ -std=gnu++11  -I"D:/PROGRA~1/R/R-42~1.1/include" -DNDEBUG  -I'D:/Program Files/R/R-4.2.1/library/cpp11/include'     -IC:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include -DHAVE_INLINE -fopenmp    -O2 -Wall  -mfpmath=sse -msse2 -mstackrealign -c dust.cpp -o dust.o

   In file included from C:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include/dust/gpu/dust_gpu.hpp:23,
                    from C:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include/dust/r/dust.hpp:15,
                    from dust.cpp:1:
   C:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include/dust/particle.hpp: In instantiation of 'dust::particle<T>::particle(dust::particle<T>::pars_type, dust::particle<T>::time_type, dust::particle<T>::rng_state_type&) [with T = odin; dust::particle<T>::pars_type = dust::pars_type<odin>; dust::particle<T>::time_type = long long unsigned int; dust::particle<T>::rng_state_type = dust::random::xoshiro_state<long long unsigned int, 4, (dust::random::scrambler)2>]':
   C:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include/dust/dust_cpu.hpp:398:36:   required from 'void dust::dust_cpu<T>::initialise(const std::vector<dust::pars_type<T> >&, dust::dust_cpu<T>::time_type, bool) [with T = odin; dust::dust_cpu<T>::time_type = long long unsigned int]'
   C:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include/dust/dust_cpu.hpp:68:5:   required from 'dust::dust_cpu<T>::dust_cpu(const std::vector<dust::pars_type<T> >&, dust::dust_cpu<T>::time_type, size_t, size_t, const std::vector<typename T::rng_state_type::int_type>&, bool, const std::vector<long long unsigned int>&) [with T = odin; dust::dust_cpu<T>::time_type = long long unsigned int; size_t = long long unsigned int; typename T::rng_state_type::int_type = long long unsigned int]'
   C:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include/dust/r/dust.hpp:48:9:   required from 'cpp11::list dust::r::dust_cpu_alloc(cpp11::list, bool, cpp11::sexp, cpp11::sexp, int, cpp11::sexp, bool, cpp11::sexp, cpp11::sexp) [with T = odin; cpp11::list = cpp11::r_vector<SEXPREC*>]'
   dust.cpp:323:64:   required from here
C:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include/dust/particle.hpp:24:26: error: no matching function for call to

   C:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include/dust/particle.hpp:24:26: error: no matching function for call to 'odin::initial(dust::particle<odin>::time_type&, dust::particle<odin>::rng_state_type&)'
        y_swap_(model_.size()) {
                             ^

   dust.cpp:31:26: note: candidate: 'std::vector<double> odin::initial(size_t)'
      std::vector<real_type> initial(size_t step) {
                             ^~~~~~~
   dust.cpp:31:26: note:   candidate expects 1 argument, 2 provided
   In file included from C:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include/dust/gpu/dust_gpu.hpp:23,
                    from C:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include/dust/r/dust.hpp:15,
                    from dust.cpp:1:
   C:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include/dust/particle.hpp: In instantiation of 'void dust::particle<T>::set_pars(dust::particle<T>::pars_type, dust::particle<T>::time_type, bool, dust::particle<T>::rng_state_type&) [with T = odin; dust::particle<T>::pars_type = dust::pars_type<odin>; dust::particle<T>::time_type = long long unsigned int; dust::particle<T>::rng_state_type = dust::random::xoshiro_state<long long unsigned int, 4, (dust::random::scrambler)2>]':
   C:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include/dust/dust_cpu.hpp:419:11:   required from 'void dust::dust_cpu<T>::initialise(const std::vector<dust::pars_type<T> >&, dust::dust_cpu<T>::time_type, bool) [with T = odin; dust::dust_cpu<T>::time_type = long long unsigned int]'
   C:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include/dust/dust_cpu.hpp:68:5:   required from 'dust::dust_cpu<T>::dust_cpu(const std::vector<dust::pars_type<T> >&, dust::dust_cpu<T>::time_type, size_t, size_t, const std::vector<typename T::rng_state_type::int_type>&, bool, const std::vector<long long unsigned int>&) [with T = odin; dust::dust_cpu<T>::time_type = long long unsigned int; size_t = long long unsigned int; typename T::rng_state_type::int_type = long long unsigned int]'
C:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include/dust/r/dust.hpp:48:9:   required from 'cpp11::list dust::r::dust_cpu_alloc(cpp11::list, bool, cpp11::sexp, cpp11::sexp, int, cpp11::sexp, bool, cpp11::sexp, cpp11

   C:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include/dust/r/dust.hpp:48:9:   required from 'cpp11::list dust::r::dust_cpu_alloc(cpp11::list, bool, cpp11::sexp, cpp11::sexp, int, cpp11::sexp, bool, cpp11::sexp, cpp11::sexp) [with T = odin; cpp11::list = cpp11::r_vector<SEXPREC*>]'
   dust.cpp:323:64:   required from here
   C:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include/dust/particle.hpp:71:18: error: passing 'const odin' as 'this' argument discards qualifiers [-fpermissive]
        if (m.size() != size()) {
   dust.cpp:28:10: note:   in call to 'size_t odin::size()'
      size_t size() {
             ^~~~
   In file included from C:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include/dust/gpu/dust_gpu.hpp:23,
                    from C:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include/dust/r/dust.hpp:15,
                    from dust.cpp:1:
   C:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include/dust/particle.hpp:74:64: error: passing 'const odin' as 'this' argument discards qualifiers [-fpermissive]
            "expected length " << size() << " but created length " <<
   dust.cpp:28:10: note:   in call to 'size_t odin::size()'
      size_t size() {
             ^~~~
   In file included from C:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include/dust/gpu/dust_gpu.hpp:23,
                    from C:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include/dust/r/dust.hpp:15,
                    from dust.cpp:1:
   C:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include/dust/particle.hpp:81:10: error: no matching function for call to 'odin::initial(dust::particle<odin>::time_type&, dust::particle<odin>::rng_state_type&)'
          y_ = model_.initial(time_, rng_state);
   dust.cpp:31:26: note: candidate: 'std::vector<double> odin::initial(size_t)'
      std::vector<real_type> initial(size_t step) {
                             ^~~~~~~
   dust.cpp:31:26: note:   candidate expects 1 argument, 2 provided

   make: *** [D:/PROGRA~1/R/R-42~1.1/etc/x64/Makeconf:229: dust.o] Error 1

   ERROR: compilation failed for package 'odin26987586'

─  removing 'C:/Users/ADMINI~1/AppData/Local/Temp/RtmpEnlFR8/devtools_install_22681092374e/odin26987586'

Error in `(function (command = NULL, args = character(), error_on_status = TRUE, …`:
! System command 'Rcmd.exe' failed
---
Exit status: 1
stdout & stderr: <printed>
---
Type .Last.error to see the more details.
richfitz commented 1 year ago

That works for me:

> sir <- odin.dust::odin_dust({
+   update(S) <- S - n_SI
+   update(I) <- I + n_SI - n_IR
+   update(R) <- R + n_IR
+   
+   n_IR <- rbinom(I, 1 - exp(-beta * I / (S + I + R)))
+   n_SI <- rbinom(S, 1 - exp(-gamma))
+   
+   initial(S) <- S_ini
+   initial(I) <- I_ini
+   initial(R) <- 0
+   
+   S_ini <- user(1000)
+   I_ini <- user(10)
+   beta <- user(0.2)
+   gamma <- user(0.1)
+ })
Loading required namespace: pkgbuild
ℹ 20 functions decorated with [[cpp11::register]]
✔ generated file cpp11.R
✔ generated file cpp11.cpp
Re-compiling odin142809c6
─  installing *source* package ‘odin142809c6’ ...
   ** using staged installation
   ** libs
   using C++ compiler: ‘g++ (Ubuntu 9.4.0-1ubuntu1~20.04.2) 9.4.0’
   g++ -std=gnu++17 -I"/usr/share/R/include" -DNDEBUG  -I'/home/rich/lib/R/library/cpp11/include' -g -Wall -Wextra -pedantic -Wmaybe-uninitialized -Wno-unused-parameter -Wno-cast-function-type -Wno-missing-field-initializers -O2   -I/home/rich/lib/R/library/dust/include -DHAVE_INLINE -fopenmp  -fpic  -g -Wall -Wextra -pedantic -Wmaybe-uninitialized -Wno-unused-parameter -Wno-cast-function-type -Wno-missing-field-initializers -O2 -c cpp11.cpp -o cpp11.o
   g++ -std=gnu++17 -I"/usr/share/R/include" -DNDEBUG  -I'/home/rich/lib/R/library/cpp11/include' -g -Wall -Wextra -pedantic -Wmaybe-uninitialized -Wno-unused-parameter -Wno-cast-function-type -Wno-missing-field-initializers -O2   -I/home/rich/lib/R/library/dust/include -DHAVE_INLINE -fopenmp  -fpic  -g -Wall -Wextra -pedantic -Wmaybe-uninitialized -Wno-unused-parameter -Wno-cast-function-type -Wno-missing-field-initializers -O2 -c dust.cpp -o dust.o
   g++ -std=gnu++17 -shared -L/usr/lib/R/lib -Wl,-Bsymbolic-functions -Wl,-z,relro -o odin142809c6.so cpp11.o dust.o -fopenmp -L/usr/lib/R/lib -lR
   installing to /tmp/RtmpbDce2V/devtools_install_77e3272c9aa5/00LOCK-file77e3745fbe27/00new/odin142809c6/libs
   ** checking absolute paths in shared objects and dynamic libraries
─  DONE (odin142809c6)
ℹ Loading odin142809c6

This bit of the error:

   C:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include/dust/particle.hpp:81:10: error: no matching function for call to 'odin::initial(dust::particle<odin>::time_type&, dust::particle<odin>::rng_state_type&)'
          y_ = model_.initial(time_, rng_state);
   dust.cpp:31:26: note: candidate: 'std::vector<double> odin::initial(size_t)'
      std::vector<real_type> initial(size_t step) {
                             ^~~~~~~
   dust.cpp:31:26: note:   candidate expects 1 argument, 2 provided

means that despite what you might think, you have an old version of dust loaded, as this was modified in https://github.com/mrc-ide/dust/pull/400

These are all errors with your system and installation, I'm afraid. Be sure to install things in a fresh session, make sure that no workspace is loaded, and check the R manuals if issues persist. Be aware that packgeVersion() reports on the installed version of a package and not the loaded version, so this can be misleading.

whzhuscu commented 1 year ago

Thanks! Now I can run the model, 

incidence <- read.csv(system.file("sir_incidence.csv", package = "mcstate"))
incidence
plot(cases ~ day, incidence,
     type = "o", xlab = "Day", ylab = "New cases", pch = 19)
history <- readRDS("F:/sir_true_history.rds")
history <- t(drop(history))
colnames(history) <- c("S", "I", "R", "cases")
history <- data.frame(t = seq_len(nrow(history)) - 1L, history)
history
matplot(history$t, history[c("S", "I", "R")], type = "l", lty = 1,
        xlab = "Day", ylab = "Variable")

sir <- odin.dust::odin_dust({
  update(S) <- S - n_SI
  update(I) <- I + n_SI - n_IR
  update(R) <- R + n_IR

  n_IR <- rbinom(I, 1 - exp(-beta * I / (S + I + R)))
  n_SI <- rbinom(S, 1 - exp(-gamma))

  initial(S) <- S_ini
  initial(I) <- I_ini
  initial(R) <- 0

  S_ini <- user(1000)
  I_ini <- user(10)
  beta <- user(0.2)
  gamma <- user(0.1)
})
dt <- 0.25
data <- mcstate::particle_filter_data(incidence, time = "day", rate = 1 / dt)

compare <- function(state, observed, pars = NULL) {
  if (is.na(observed$cases)) {
    return(NULL)
  }
  exp_noise <- 1e6
  incidence_modelled <- state[1, , drop = TRUE]
  incidence_observed <- observed$cases
  lambda <- incidence_modelled +
    rexp(n = length(incidence_modelled), rate = exp_noise)
  dpois(x = incidence_observed, lambda = lambda, log = TRUE)
}
index <- function(info) {
  list(run = 5L, state = 1:3)
}

proposal_kernel <- rbind(c(0.00057, 0.00034), c(0.00034, 0.00026))
pars <- mcstate::pmcmc_parameters$new(
  list(mcstate::pmcmc_parameter("beta", 0.2, min = 0, max = 1,
                                prior = function(p) log(1e-10)),
       mcstate::pmcmc_parameter("gamma", 0.1, min = 0, max = 1,
                                prior = function(p) log(1e-10))),
  proposal = proposal_kernel)
pars
n_steps <- 500
n_particles <- 100
n_threads <- dust::dust_openmp_threads()

p1 <- mcstate::particle_filter$new(data, sir, n_particles, compare,
                                   index = index,
                                   n_threads = n_threads)
control1 <- mcstate::pmcmc_control(n_steps, save_trajectories = TRUE,
                                   save_state = TRUE, save_restart = 40)
res1 <- mcstate::pmcmc(pars, p1, control = control1)

I also modify the state:

index <- function(info) {
  list(run = 5L, state = 1:3)
}

But there is error, could you please give me some suggestion? Thank you.

> res1 <- mcstate::pmcmc(pars, p1, control = control1)
Error: All elements of 'index' must lie in [1, 3]
whzhuscu commented 1 year ago

An example in your paper is,

sir <- dust::dust_example("sir")

but I do not konw what the model looks like. So it is difficult for me to rewrite the code.

> sir
<dust> object generator
  Public:
    initialize: function (pars, time, n_particles, n_threads = 1L, seed = NULL, 
    name: function () 
    param: function () 
    run: function (time_end) 
    simulate: function (time_end) 
    set_index: function (index) 
    index: function () 
    ode_control: function () 
    ode_statistics: function () 
    n_threads: function () 
    n_state: function () 
    n_particles: function () 
    n_particles_each: function () 
    shape: function () 
    update_state: function (pars = NULL, state = NULL, time = NULL, set_initial_state = NULL, 
    state: function (index = NULL) 
    time: function () 
    set_stochastic_schedule: function (time) 
    reorder: function (index) 
    resample: function (weights) 
    info: function () 
    pars: function () 
    rng_state: function (first_only = FALSE, last_only = FALSE) 
    set_rng_state: function (rng_state) 
    has_openmp: function () 
    has_gpu_support: function (fake_gpu = FALSE) 
    has_compare: function () 
    real_size: function () 
    time_type: function () 
    rng_algorithm: function () 
    uses_gpu: function (fake_gpu = FALSE) 
    n_pars: function () 
    set_n_threads: function (n_threads) 
    set_data: function (data, shared = FALSE) 
    compare_data: function () 
    filter: function (time_end = NULL, save_trajectories = FALSE, time_snapshot = NULL, 
    gpu_info: function () 
  Private:
    pars_: NULL
    pars_multi_: NULL
    index_: NULL
    info_: NULL
    n_threads_: NULL
    n_particles_: NULL
    n_particles_each_: NULL
    shape_: NULL
    ptr_: NULL
    gpu_config_: NULL
    ode_control_: NULL
    methods_: NULL
    param_: list
    reload_: NULL
  Parent env: <environment: namespace:dust>
  Locked objects: TRUE
  Locked class: FALSE
  Portable: TRUE
richfitz commented 1 year ago

Hi,

I am afraid that issues are not the right place to seek technical support with the packages, and in general - like many open source packages - we cannot provide technical support directly. Issues should be used for potential problems with the package.

You have opened issues on several different repos, and repeatedly opened and closed them which makes finding the right place to respond to very hard.

The error here:

> res1 <- mcstate::pmcmc(pars, p1, control = control1)
Error: All elements of 'index' must lie in [1, 3]

I explained on one of your other issues - you are trying to access elements 1:5 from the model but your model only has 3 variables, which is what the error is telling you. Please see the documentation for mcstate and dust where this is discussed.

but I do not konw what the model looks like. So it is difficult for me to rewrite the code.

Please see the dust documentation, where this is discussed: https://mrc-ide.github.io/dust/articles/dust.html

More documentation around the set of packages can be found here: https://mrc-ide.github.io/odin/articles/guide.html

Perhaps try working through some of the examples in the mcstate or odin.dust documentation. These are all run automatically as part of our testing, so they do work - if you run into issues please see the R documentation.

whzhuscu commented 1 year ago

Thanks. I am sorry, my Internet connection is so bad that I have repeatedly opened and closed issues.

richfitz commented 1 year ago

No worries, I just get lots of emails from github when it happens and it looks like someone trying to get my attention.

These are the specific docs that you probably want to work through, I suggest that these will be easier than the paper code:

whzhuscu commented 1 year ago

Thanks for your advice, I have been learning your and your team's code. But this example of using PMCMC to estimate parameters of stochastic SIR model is a model read directly with this statement.

sir <- dust::dust_example("sir")

I don't know the specific form of code of this model, and I don't know how to use '_odin.dust::odindust' to build a stochastic SEIR model. I wonder if it is convenient for you to give me some suggestions. Sorry, I am a little poor at coding. Thank you very much! And the code has the error:

res1 <- mcstate::pmcmc(pars, p1, control = control1)
Error: All elements of 'index' must lie in [1, 3] 
14. private$methods_$set_index(private$ptr_, index) at dust.R#95
13. model$set_index(index_data$run) 
12. initialize(...) 
11. particle_filter_state$new(pars, self$model, private$last_model[[1]], 
    private$data, private$data_split, private$times, self$n_particles, 
    self$has_multiple_parameters, private$n_threads, private$initial, 
    private$index, private$compare, private$constant_log_likelihood,  ... 
10. self$run_begin(pars, save_history, save_restart, min_log_likelihood = min_log_likelihood) 
9. filter_run_simple(self, private, pars, save_history, save_restart, 
    min_log_likelihood) 
8. filter_run(self, private, pars, save_history, save_restart, min_log_likelihood) 
7. private$filter$run(private$pars$model(p), private$control$save_trajectories, 
    private$control$save_restart, min_log_likelihood) 
6. private$run_filter(private$curr_pars) 
5. initialize(...) 
4. pmcmc_state$new(pars, initial, filter, control) 
3. pmcmc_run_chain(pars, initial[[i]], filter, control, NULL) 
2. pmcmc_multiple_series(pars, initial, filter, control) 
1. mcstate::pmcmc(pars, p1, control = control1)