mjg211 / phaseR

Development version of phaseR, an R package for phase plane analysis of one- and two-dimensional autonomous ODE systems
https://doi.org/10.32614/RJ-2014-023
Other
15 stars 3 forks source link

Error when writing equations in matrix form #14

Open juanrocha opened 9 months ago

juanrocha commented 9 months ago

Thanks for developing and maintaining the package!

I run into the following error when trying to plot a phase diagram on a 2-dim system. The function is:

## Pollution model:
pollution <- function(t, y, params){
  with(as.list(c(y, params)), {
    x = numeric(length = n)
    x = y
    pollutant <- u - s*x + v * (x^alpha/(z^alpha + x^alpha))
    outflow <-  (A_ij * delta_ij) %*% x
    inflow <-  t(A_ij * delta_ij) %*% x

    dx <- pollutant + (inflow - outflow)

    return(list(c(dx)))
  })
}

With the following parameters:

n <- 1 # number of systems
delta_ij <- matrix(runif(n^2, min = 0.02, max = 0.05), ncol = n)
A_ij <- matrix(rbinom(n^2, 1, prob = 0.3), n, n)

## Parameters: 
params <- list(
    n = n,                             # number of systems (e.g. lakes)
    u = rep(0.5, n),             # pollution load from humans
    s = rep(2.5, n),             # internal loss rate (sedimentation)
    v = rep(10, n),              # max level of internal nutrient release
    z = rep(2.2, n),             # threshold
    alpha = 4 ,                    # sharpness of the shift
    delta_ij = delta_ij,        # matrix of difussion terms
    A_ij = A_ij                     # adjacency matrix
)

phaseR works on the one-dimensional case, eg:

fld <- flowField(pollution, xlim = c(0,5), ylim = c(-2,7), parameters = params,
          system = "one.dim", add = FALSE)
nll <- nullclines(pollution, xlim = c(0,5), ylim = c(-2,7), parameters = params,
          system = "one.dim", add = TRUE)
trj <- trajectory(pollution, y0 = c(-0.5, 1, 3, 5), tlim = c(0,5), parameters = params,
          system = "one.dim", add = TRUE)

But as soon as I change n = 2 I get the following error:

fld <- flowField(pollution, xlim = c(0,5), ylim = c(-2,7), system = "two.dim", add = FALSE, parameters = params)

Error in (A_ij * delta_ij) %*% x : non-conformable arguments

I don't quite understand the origin of the error. The system works with n = 1, 2, or 100s when doing numerical simulations with deSolve, but for some reason fails on the two-dimensional case in phaseR.

If I rewrite the function in a non-matrix form such as:

pollution2 <- function(t,y,params){
    with(as.list(c(y, params)), {
    x = numeric(length = n)
    x = y
    dx <- u[1] - s[1]*x[1] + v[1] * (x[1]^alpha/(z[1]^alpha + x[1]^alpha)) + (A_ij[2,1] * delta_ij[2,1]) %*% x[1] - (A_ij[1,2] * delta_ij[1,2]) %*% x[2]

    dy <- u[2] - s[2]*x[2] + v[2] * (x[2]^alpha/(z[2]^alpha + x[2]^alpha)) + (A_ij[1,2] * delta_ij[1,2]) %*% x[2] - (A_ij[1,2] * delta_ij[1,2]) %*% x[1]

    return(list(c(dx, dy)))
  })
}

Trying to get the flow field throughs a different error:

fld <- flowField(pollution2, xlim = c(0,5), ylim = c(-2,7), system = "two.dim", add = FALSE, parameters = params)

Error in if (any(dx[i, j] != 0, dy[i, j] != 0)) { : missing value where TRUE/FALSE needed

Which seems similar to issue #13 Any help is greatly appreciated!

Below my session info if useful:

session_info()
─ Session info ─────────────────────────────────────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.2.2 (2022-10-31)
 os       macOS Ventura 13.5.2
 system   x86_64, darwin17.0
 ui       RStudio
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       Europe/Stockholm
 date     2023-09-23
 rstudio  2023.06.0+421 Mountain Hydrangea (desktop)
 pandoc   2.14.2 @ /usr/local/bin/pandoc

─ Packages ─────────────────────────────────────────────────────────────────────────────────────────────────
 package        * version date (UTC) lib source
 assertthat       0.2.1   2019-03-21 [1] CRAN (R 4.2.0)
 backports        1.4.1   2021-12-13 [1] CRAN (R 4.2.0)
 broom            1.0.2   2022-12-15 [1] CRAN (R 4.2.2)
 cachem           1.0.7   2023-02-24 [1] CRAN (R 4.2.0)
 callr            3.7.3   2022-11-02 [1] CRAN (R 4.2.0)
 cellranger       1.1.0   2016-07-27 [1] CRAN (R 4.2.0)
 cli              3.6.1   2023-03-23 [1] CRAN (R 4.2.0)
 coda             0.19-4  2020-09-30 [1] CRAN (R 4.2.0)
 colorspace       2.1-0   2023-01-23 [1] CRAN (R 4.2.0)
 crayon           1.5.2   2022-09-29 [1] CRAN (R 4.2.0)
 DBI              1.1.3   2022-06-18 [1] CRAN (R 4.2.0)
 dbplyr           2.2.1   2022-06-27 [1] CRAN (R 4.2.0)
 deSolve        * 1.34    2022-10-22 [1] CRAN (R 4.2.0)
 devtools       * 2.4.5   2022-10-11 [1] CRAN (R 4.2.0)
 digest           0.6.33  2023-07-07 [1] CRAN (R 4.2.0)
 dplyr          * 1.1.1   2023-03-22 [1] CRAN (R 4.2.0)
 ellipsis         0.3.2   2021-04-29 [1] CRAN (R 4.2.0)
 fansi            1.0.4   2023-01-22 [1] CRAN (R 4.2.0)
 fastmap          1.1.1   2023-02-24 [1] CRAN (R 4.2.0)
 forcats        * 0.5.2   2022-08-19 [1] CRAN (R 4.2.0)
 fs               1.6.3   2023-07-20 [1] CRAN (R 4.2.0)
 gargle           1.2.1   2022-09-08 [1] CRAN (R 4.2.0)
 generics         0.1.3   2022-07-05 [1] CRAN (R 4.2.0)
 ggplot2        * 3.4.2   2023-04-03 [1] CRAN (R 4.2.0)
 glue             1.6.2   2022-02-24 [1] CRAN (R 4.2.0)
 googledrive      2.0.0   2021-07-08 [1] CRAN (R 4.2.0)
 googlesheets4    1.0.1   2022-08-13 [1] CRAN (R 4.2.0)
 gtable           0.3.3   2023-03-21 [1] CRAN (R 4.2.0)
 haven            2.5.1   2022-08-22 [1] CRAN (R 4.2.0)
 hms              1.1.2   2022-08-19 [1] CRAN (R 4.2.0)
 htmltools        0.5.5   2023-03-23 [1] CRAN (R 4.2.0)
 htmlwidgets      1.6.0   2022-12-15 [1] CRAN (R 4.2.2)
 httpuv           1.6.9   2023-02-14 [1] CRAN (R 4.2.0)
 httr             1.4.5   2023-02-24 [1] CRAN (R 4.2.0)
 igraph           1.3.5   2022-09-22 [1] CRAN (R 4.2.0)
 jsonlite         1.8.7   2023-06-29 [1] CRAN (R 4.2.0)
 knitr            1.43    2023-05-25 [1] CRAN (R 4.2.0)
 later            1.3.0   2021-08-18 [1] CRAN (R 4.2.0)
 lattice          0.20-45 2021-09-22 [1] CRAN (R 4.2.2)
 lifecycle        1.0.3   2022-10-07 [1] CRAN (R 4.2.0)
 lubridate        1.9.2   2023-02-10 [1] CRAN (R 4.2.0)
 magrittr         2.0.3   2022-03-30 [1] CRAN (R 4.2.0)
 memoise          2.0.1   2021-11-26 [1] CRAN (R 4.2.0)
 mime             0.12    2021-09-28 [1] CRAN (R 4.2.0)
 miniUI           0.1.1.1 2018-05-18 [1] CRAN (R 4.2.0)
 modelr           0.1.10  2022-11-11 [1] CRAN (R 4.2.0)
 munsell          0.5.0   2018-06-12 [1] CRAN (R 4.2.0)
 network          1.18.0  2022-10-06 [1] CRAN (R 4.2.0)
 phaseR         * 2.2.1   2022-09-02 [1] CRAN (R 4.2.0)
 pillar           1.9.0   2023-03-22 [1] CRAN (R 4.2.0)
 pkgbuild         1.4.0   2022-11-27 [1] CRAN (R 4.2.0)
 pkgconfig        2.0.3   2019-09-22 [1] CRAN (R 4.2.0)
 pkgload          1.3.2.1 2023-07-08 [1] CRAN (R 4.2.0)
 prettyunits      1.1.1   2020-01-24 [1] CRAN (R 4.2.0)
 processx         3.8.2   2023-06-30 [1] CRAN (R 4.2.0)
 profvis          0.3.7   2020-11-02 [1] CRAN (R 4.2.0)
 promises         1.2.0.1 2021-02-11 [1] CRAN (R 4.2.0)
 ps               1.7.5   2023-04-18 [1] CRAN (R 4.2.0)
 purrr          * 1.0.2   2023-08-10 [1] CRAN (R 4.2.0)
 R6               2.5.1   2021-08-19 [1] CRAN (R 4.2.0)
 Rcpp             1.0.10  2023-01-22 [1] CRAN (R 4.2.0)
 readr          * 2.1.3   2022-10-01 [1] CRAN (R 4.2.0)
 readxl           1.4.1   2022-08-17 [1] CRAN (R 4.2.0)
 remotes          2.4.2.1 2023-07-18 [1] CRAN (R 4.2.0)
 reprex           2.0.2   2022-08-17 [1] CRAN (R 4.2.0)
 rlang            1.1.1   2023-04-28 [1] CRAN (R 4.2.0)
 rstudioapi       0.14    2022-08-22 [1] CRAN (R 4.2.0)
 rvest            1.0.3   2022-08-19 [1] CRAN (R 4.2.0)
 scales           1.2.1   2022-08-20 [1] CRAN (R 4.2.0)
 sessioninfo      1.2.2   2021-12-06 [1] CRAN (R 4.2.0)
 shiny            1.7.4   2022-12-15 [1] CRAN (R 4.2.2)
 statnet.common   4.7.0   2022-09-08 [1] CRAN (R 4.2.0)
 stringi          1.7.12  2023-01-11 [1] CRAN (R 4.2.0)
 stringr        * 1.5.0   2022-12-02 [1] CRAN (R 4.2.0)
 tibble         * 3.2.1   2023-03-20 [1] CRAN (R 4.2.0)
 tidyr          * 1.3.0   2023-01-24 [1] CRAN (R 4.2.0)
 tidyselect       1.2.0   2022-10-10 [1] CRAN (R 4.2.0)
 tidyverse      * 1.3.2   2022-07-18 [1] CRAN (R 4.2.0)
 timechange       0.2.0   2023-01-11 [1] CRAN (R 4.2.0)
 tzdb             0.3.0   2022-03-28 [1] CRAN (R 4.2.0)
 urlchecker       1.0.1   2021-11-30 [1] CRAN (R 4.2.0)
 usethis        * 2.1.6   2022-05-25 [1] CRAN (R 4.2.0)
 utf8             1.2.3   2023-01-31 [1] CRAN (R 4.2.0)
 vctrs            0.6.3   2023-06-14 [1] CRAN (R 4.2.0)
 waydown        * 1.1.0   2021-03-17 [1] CRAN (R 4.2.0)
 withr            2.5.0   2022-03-03 [1] CRAN (R 4.2.0)
 xfun             0.40    2023-08-09 [1] CRAN (R 4.2.0)
 xml2             1.3.5   2023-07-06 [1] CRAN (R 4.2.0)
 xtable           1.8-4   2019-04-21 [1] CRAN (R 4.2.0)

 [1] /Library/Frameworks/R.framework/Versions/4.2/Resources/library