Code Monkey home page Code Monkey logo

phaser's People

Contributors

burgerga avatar mjg211 avatar tomicapretto avatar

Stargazers

avatar avatar avatar avatar avatar avatar avatar avatar avatar avatar avatar avatar avatar avatar avatar

Watchers

avatar avatar

Forkers

burgerga tomicapretto bbolker

phaser's Issues

Method "ode45" not adapted for some models trajectories

phaseR/R/trajectory.R

Line 209 in 4f58357

method <- ifelse(tstep > 0, "ode45", "lsoda")

First of thank you for your time-saving package. I give a reproducible exemple below: for some initial conditions on this model, method "ode45" computing time for the trajectory is exponentially longer than method "euler". I am not sure whether the choice of "ode45" was guided by some particular reason? Maybe you could add the possibility to provide a method in the call to trajectory?

library(tictoc)
library(deSolve)

params <- c(s = 0.35, gamma = 2.5, n = 1.5, alpha = 3, beta = 2)

model <- function(t, y, parameters){
  with(as.list(parameters),{
    g <- y[1] 
    gw <- y[2]
    dg <- gamma*s*(1/gw-1/g)
    dgw <- -gw^2/s*(alpha*s*(1/g-1/gw) - beta*(g-n))
    
    list(c(dg,dgw))
  })
}
yini <- c(g = 0.27, gw = 1)
times <- seq(0, 5, by = 0.1)
tic()
out <- ode(yini, times, model, params, method = "ode45")
toc()
24.501 sec elapsed

tic()
out <- ode(yini, times, model, params, method = "euler")
toc()
0.003 sec elapsed

Possible typo in RJournal manuscript

Hi! Thanks for this fantastic package.

Not sure if this is the place for this, but I think there might be a typo in the associated R Journal manuscript (found here)

image

I think this is the other way around, right? (e.g. dx/dt = -x is stable at x=0, and k<0; dx/dt = x is unstable at x=0, and k>0)

Error when writing equations in matrix form

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

Missing value where TRUE/FALSE needed

I am trying to draw the phase diagram of the Solow economic growth model. Using y^a with a that is a non-natural number, I get:

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

The code I use is:

solow <- function(t, y, parameters) {
  list(y^0.5)
}
example2_flowField     <- flowField(solow,
                                    xlim   = c(0, 4),
                                    ylim   = c(0, 4),
                                    system = "one.dim",
                                    add    = FALSE,
                                    xlab   = "t")

Can somebody help me?

deriv function has inflexible `parameters` argument

Sorry for the vague title, it is a mistake in my own code:

library(phaseR)
#> -------------------------------------------------------------------------------
#> phaseR: Phase plane analysis of one- and two-dimensional autonomous ODE systems
#> -------------------------------------------------------------------------------
#> 
#> v.2.1: For an overview of the package's functionality enter: ?phaseR
#> 
#> For news on the latest updates enter: news(package = "phaseR")
ex2.alt <- function(t, state, parameters) {
    with(as.list(c(state, parameters)), {
      dA = A*(1 - A)*(2 - A)
      list(c(dA))
    })
}
example2_stability_1 <- stability(ex2.alt, ystar = 0, system = "one.dim", 
                                  state.names = "A")
#> discriminant = 2, classification = Unstable
                                  
ex2.alt.alt <- function(t, state, pars) {
    with(as.list(c(state, pars)), {
      dA = A*(1 - A)*(2 - A)
      list(c(dA))
    })
}
example2_stability_1 <- stability(ex2.alt.alt, ystar = 0, system = "one.dim", 
                                  state.names = "A")
#> Error in deriv(0, stats::setNames(ystar + h, state.names), parameters = parameters): unused argument (parameters = parameters)

Created on 2019-09-08 by the reprex package (v0.3.0)

Apparently, it is not possible to change the name of the parameters argument in your deriv
The offending lines are https://github.com/mjg211/phaseR/blob/615b31c7ac/R/stability.R#L100-L101
Solution is easy: just change parameters = parameters into parameters (will make a pull request when I have time, but if you're earlier maybe you can take care of this when resolving #6 )

vignette("phaseR") does not work

Just a (self?) reminder that vignette("phaseR"), as suggested in ?phaseR, does not work.
The call that work is vignette("my-vignette", package = "phaseR")

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.