From a3eaa31f476aa5d86621ec2eb7242642da3be890 Mon Sep 17 00:00:00 2001 From: Neves-P Date: Tue, 7 Mar 2023 13:45:30 +0100 Subject: [PATCH 01/13] Cpp14 Makevars --- DESCRIPTION | 7 +++---- src/Makevars | 2 ++ src/Makevars.win | 2 ++ src/config.h | 7 +++++-- 4 files changed, 12 insertions(+), 6 deletions(-) create mode 100644 src/Makevars create mode 100644 src/Makevars.win diff --git a/DESCRIPTION b/DESCRIPTION index fb8b8f83..4fcee504 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: DAISIE Type: Package Title: Dynamical Assembly of Islands by Speciation, Immigration and Extinction -Version: 4.3.1.9000 -Date: 2023-02-10 +Version: 4.3.2 +Date: 2023-02-24 Depends: R (>= 4.1.0) biocViews: Imports: @@ -110,8 +110,7 @@ Description: Simulates and computes the (maximum) likelihood of a dynamical model of island biota assembly through speciation, immigration and extinction. See e.g. Valente et al. 2015. Ecology Letters 18: 844-852, . -NeedsCompilation: yes -SystemRequirements: C++14 +NeedsCompilation: yesc Encoding: UTF-8 VignetteBuilder: knitr URL: https://github.com/rsetienne/DAISIE diff --git a/src/Makevars b/src/Makevars new file mode 100644 index 00000000..1fada8d1 --- /dev/null +++ b/src/Makevars @@ -0,0 +1,2 @@ +CXX_STD = CXX14 +PKG_CPPFLAGS = -D_HAS_AUTO_PTR_ETC=0 diff --git a/src/Makevars.win b/src/Makevars.win new file mode 100644 index 00000000..1fada8d1 --- /dev/null +++ b/src/Makevars.win @@ -0,0 +1,2 @@ +CXX_STD = CXX14 +PKG_CPPFLAGS = -D_HAS_AUTO_PTR_ETC=0 diff --git a/src/config.h b/src/config.h index 9dea9a57..a14c60c8 100644 --- a/src/config.h +++ b/src/config.h @@ -1,16 +1,19 @@ -// [[Rcpp::plugins(cpp14)]] - #ifndef ODEINT_CONFIG_H_INCLUDED #define ODEINT_CONFIG_H_INCLUDED +// [[Rcpp::plugins(cpp14)]] + // Special case to make use of some steppers that would include // boost/functional.hpp +// moved to Makevars[.win] +/* #if __cplusplus >= 201703L #ifdef _HAS_AUTO_PTR_ETC #undef _HAS_AUTO_PTR_ETC #endif #define _HAS_AUTO_PTR_ETC 0 #endif + */ // Special case to make use of some steppers that would include // boost/get_pointer.hpp From ca2c7c295b3c6add38f545c3998fd5e7606b97c4 Mon Sep 17 00:00:00 2001 From: Neves-P Date: Tue, 7 Mar 2023 13:53:30 +0100 Subject: [PATCH 02/13] Increment version number to 4.3.3 [run ci] --- DESCRIPTION | 2 +- NEWS.md | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 4fcee504..1ab691c1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: DAISIE Type: Package Title: Dynamical Assembly of Islands by Speciation, Immigration and Extinction -Version: 4.3.2 +Version: 4.3.3 Date: 2023-02-24 Depends: R (>= 4.1.0) biocViews: diff --git a/NEWS.md b/NEWS.md index 5aff1939..6effaa8a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,4 @@ -# DAISIE (development version) +# DAISIE 4.3.3 # DAISIE 4.3.1 From 92d86481566f8f18e91e1cba9e99f525351bf717 Mon Sep 17 00:00:00 2001 From: Hanno Hildenbrandt Date: Thu, 9 Mar 2023 11:07:20 +0100 Subject: [PATCH 03/13] odeint valgrind issue fixed --- .Rbuildignore | 3 + .vscode/c_cpp_properties.json | 47 +++ .vscode/genenv.R | 19 + .vscode/launch.json | 40 ++ .vscode/tasks.json | 13 + NAMESPACE | 1 - source_odeint.R | 242 +++++++++++ src/DAISIE_loglik_rhs_FORTRAN.f95 | 248 +++++------ src/DAISIE_odeint.h | 1 + src/config.h | 17 +- src/patched_bulrisch_stoer.hpp | 661 ++++++++++++++++++++++++++++++ 11 files changed, 1160 insertions(+), 132 deletions(-) create mode 100644 .vscode/c_cpp_properties.json create mode 100644 .vscode/genenv.R create mode 100644 .vscode/launch.json create mode 100644 .vscode/tasks.json create mode 100644 source_odeint.R create mode 100644 src/patched_bulrisch_stoer.hpp diff --git a/.Rbuildignore b/.Rbuildignore index 5d79eddc..5954eef4 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -16,3 +16,6 @@ ^\.covrignore$ ^LICENSE$ ^\.zenodo\.json$ + +^\.vscode$ +^source_odeint\.R diff --git a/.vscode/c_cpp_properties.json b/.vscode/c_cpp_properties.json new file mode 100644 index 00000000..19fc441b --- /dev/null +++ b/.vscode/c_cpp_properties.json @@ -0,0 +1,47 @@ +// c_cpp_properties.json +// +// C/C++ extension config file +// Some settings to augment IntelliSense for C/C++ files +{ + "configurations": [ + { + "name": "Linux gdb", + "includePath": [ + // recurse into workspace, picks up ./src/ and ./inst/include/ + "${workspaceFolder}/**", + // R headers + "${env:HOME}/opt/bin/Rroot/lib/R/include", + // includes used by our package + "${env:HOME}/opt/bin/Rlibrary/Rcpp/include" + // boost: "${env:HOME}/opt/bin/Rlibrary/BH/include" + ], + "defines": [], + // absolute path to compiler + "compilerPath": "/usr/bin/gcc", + // "compilerPath": "/usr/bin/clang" + "cStandard": "c11", + // should match [[Rcpp::plugins(cpp17)]] + "cppStandard": "c++17" + }, + { + "name": "Linux lldb", + "includePath": [ + // recurse into workspace, picks up ./src/ and ./inst/include/ + "${workspaceFolder}/**", + // R headers + "${env:HOME}/opt/bin/Rroot/lib/R/include", + // includes used by our package + "${env:HOME}/opt/bin/Rlibrary/Rcpp/include" + // boost: "${env:HOME}/opt/bin/Rlibrary/BH/include" + ], + "defines": [], + // absolute path to compiler + "compilerPath": "/usr/bin/clang", + // "compilerPath": "/usr/bin/clang" + "cStandard": "c11", + // should match [[Rcpp::plugins(cpp17)]] + "cppStandard": "c++17" + } + ], + "version": 4 +} diff --git a/.vscode/genenv.R b/.vscode/genenv.R new file mode 100644 index 00000000..2cd0a945 --- /dev/null +++ b/.vscode/genenv.R @@ -0,0 +1,19 @@ +# genenv.R +# +# Colects environment variables from R, required by the R *binary*. +# They are normally set by the R *script* which is circumvented if +# we need to run the binary under the debugger as we want to debug R +# not bash... +# +# This script is configured as a task in 'tasks.json' and referenced within +# 'launch.json' as the 'preLaunchTask'. + +env <- Sys.getenv() +envnames <- names(env) +rnames <- envnames[startsWith(envnames, "R_")] +cached_names <- rnames +ld_lib_path <- Sys.getenv("LD_LIBRARY_PATH") +if (ld_lib_path != "") { + cached_names <- c("LD_LIBRARY_PATH", rnames) +} +writeLines(paste0(cached_names, "=", env[cached_names]), ".vscode/.env") diff --git a/.vscode/launch.json b/.vscode/launch.json new file mode 100644 index 00000000..665df8ae --- /dev/null +++ b/.vscode/launch.json @@ -0,0 +1,40 @@ +// launch.json +// +// launch the R *binary* under gdb and run devtools::test() +{ + // Use IntelliSense to learn about possible attributes. + // Hover to view descriptions of existing attributes. + // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387 + "version": "0.2.0", + "configurations": [ + { + "name": "(gbd) devtools::test()", + "type": "cppdbg", + "request": "launch", + // The binary, not the script + "program": "${env:HOME}/opt/bin/Rroot/lib/R/bin/exec/R", + "args": [ + "--vanilla", + "-e", + "devtools::test()" + ], + "stopAtEntry": false, + // needs to be generated, see below + "envFile": "${workspaceFolder}/.vscode/.env", + "cwd": "${workspaceFolder}", + "externalConsole": false, + "MIMode": "gdb", + "setupCommands": [ + { + "description": "Enable pretty-printing for gdb", + "text": "-enable-pretty-printing", + "ignoreFailures": true + } + ], + // 'R' is a script that sets a ton of environment variables + // required by the R binary. This task emulates that part of + // the R script: + "preLaunchTask": "genenv" + } + ] +} diff --git a/.vscode/tasks.json b/.vscode/tasks.json new file mode 100644 index 00000000..e11a2f20 --- /dev/null +++ b/.vscode/tasks.json @@ -0,0 +1,13 @@ +{ + // See https://go.microsoft.com/fwlink/?LinkId=733558 + // for the documentation about the tasks.json format + "version": "2.0.0", + "tasks": [ + { + "label": "genenv", + "type": "shell", + "command": "Rscript ${workspaceFolder}/.vscode/genenv.R", + "problemMatcher": [] + } + ] +} diff --git a/NAMESPACE b/NAMESPACE index 3cf5beee..9f793cfe 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -42,7 +42,6 @@ export(create_pars) export(create_trait_pars) export(create_trait_pars_2K) export(daisie_odeint_cs) -export(daisie_odeint_iw) import(Rcpp) importFrom(doParallel,registerDoParallel) importFrom(foreach,"%dopar%") diff --git a/source_odeint.R b/source_odeint.R new file mode 100644 index 00000000..8a6ea202 --- /dev/null +++ b/source_odeint.R @@ -0,0 +1,242 @@ +# source_odeint.R +# +# partial odeint test suite from tests/testthat/test_odeint.R +# Run with: +# +# R -d "valgrind --tool=memcheck --leak-check=full --error-exitcode=1" --vanilla -e "source('source_odeint.R')" &> source_odeint.log + + +library(DAISIE) + +utils::data(Galapagos_datalist_2types) +pars1 <- c( +0.195442017, +0.087959583, +Inf, +0.002247364, +0.873605049, +3755.202241, +8.909285094, +14.99999923, +0.002247364, +0.873605049, +0.163 +) +pars2 <- c(40, 11, 0, 1) +methode <- 'deSolve_R::lsodes' +loglik_lsodes_R <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) +methode <- 'deSolve_R::lsoda' +loglik_lsoda_R <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) +methode <- 'lsodes' +loglik_lsodes <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) +methode <- 'lsoda' +loglik_lsoda <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) +methode <- 'odeint::runge_kutta_cash_karp54' +loglik_rkck54 <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) +methode <- 'odeint::runge_kutta_fehlberg78' +loglik_rkf78 <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) +methode <- 'odeint::runge_kutta_dopri5' +loglik_rkd5 <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) +methode <- 'odeint::bulirsch_stoer' +loglik_bs <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) +methode <- 'odeint::rosenbrock4' +loglik_rb <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) +methode <- 'odeint::adams_bashforth_moulton_1' +DAISIE_CS_max_steps(100000000) +DAISIE_abm_factor(0.000001) +loglik_abm <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) + + +==2373== Memcheck, a memory error detector +==2373== Copyright (C) 2002-2022, and GNU GPL'd, by Julian Seward et al. +==2373== Using Valgrind-3.19.0 and LibVEX; rerun with -h for copyright info +==2373== Command: /home/hanno/opt/Rval/lib/R/bin/exec/R --vanilla -e source('source_odeint.R') +==2373== + +R Under development (unstable) (2023-03-07 r83950) -- "Unsuffered Consequences" +Copyright (C) 2023 The R Foundation for Statistical Computing +Platform: x86_64-pc-linux-gnu (64-bit) + +R is free software and comes with ABSOLUTELY NO WARRANTY. +You are welcome to redistribute it under certain conditions. +Type 'license()' or 'licence()' for distribution details. + +R is a collaborative project with many contributors. +Type 'contributors()' for more information and +'citation()' on how to cite R or R packages in publications. + +Type 'demo()' for some demos, 'help()' for on-line help, or +'help.start()' for an HTML browser interface to help. +Type 'q()' to quit R. + +> source('source_odeint.R') +Status of colonist: 0, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -0.007835 +Status of colonist: 0, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: -0.008967 +Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 +Status of colonist: 4, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.497364 +Status of colonist: 2, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: 1.545883 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -13.576394 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.966530 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -8.752420 +Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.981045 +Parameters: +lambda^c, mu, K, gamma, lambda^a +0.195442, 0.087960, Inf, 0.002247, 0.873605 +Loglikelihood: -61.702819 + +Status of colonist: 0, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -0.007835 +Status of colonist: 0, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: -0.008967 +Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 +Status of colonist: 4, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.497364 +Status of colonist: 2, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: 1.545883 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -13.576394 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.966530 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -8.752420 +Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.981045 +Parameters: +lambda^c, mu, K, gamma, lambda^a +0.195442, 0.087960, Inf, 0.002247, 0.873605 +Loglikelihood: -61.702819 + +Status of colonist: 0, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -0.007835 +Status of colonist: 0, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: -0.008967 +Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 +Status of colonist: 4, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.497364 +Status of colonist: 2, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: 1.545883 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -13.576394 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.966530 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -8.752420 +Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.981045 +Parameters: +lambda^c, mu, K, gamma, lambda^a +0.195442, 0.087960, Inf, 0.002247, 0.873605 +Loglikelihood: -61.702819 + +Status of colonist: 0, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -0.007835 +Status of colonist: 0, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: -0.008967 +Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 +Status of colonist: 4, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.497364 +Status of colonist: 2, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: 1.545883 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -13.576394 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.966530 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -8.752420 +Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.981045 +Parameters: +lambda^c, mu, K, gamma, lambda^a +0.195442, 0.087960, Inf, 0.002247, 0.873605 +Loglikelihood: -61.702819 + +Status of colonist: 0, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -0.007835 +Status of colonist: 0, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: -0.008967 +Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 +Status of colonist: 4, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.497364 +Status of colonist: 2, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: 1.545883 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -13.576394 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.966530 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -8.752420 +Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.981045 +Parameters: +lambda^c, mu, K, gamma, lambda^a +0.195442, 0.087960, Inf, 0.002247, 0.873605 +Loglikelihood: -61.702819 + +Status of colonist: 0, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -0.007835 +Status of colonist: 0, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: -0.008967 +Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 +Status of colonist: 4, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.497364 +Status of colonist: 2, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: 1.545883 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -13.576394 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.966530 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -8.752420 +Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.981045 +Parameters: +lambda^c, mu, K, gamma, lambda^a +0.195442, 0.087960, Inf, 0.002247, 0.873605 +Loglikelihood: -61.702819 + +Status of colonist: 0, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -0.007835 +Status of colonist: 0, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: -0.008967 +Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 +Status of colonist: 4, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.497364 +Status of colonist: 2, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: 1.545883 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -13.576394 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.966530 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -8.752420 +Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.981045 +Parameters: +lambda^c, mu, K, gamma, lambda^a +0.195442, 0.087960, Inf, 0.002247, 0.873605 +Loglikelihood: -61.702819 + +Status of colonist: 0, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -0.007835 +Status of colonist: 0, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: -0.008967 +Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 +Status of colonist: 4, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.497364 +Status of colonist: 2, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: 1.545883 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -13.576394 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.966530 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -8.752420 +Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.981045 +Parameters: +lambda^c, mu, K, gamma, lambda^a +0.195442, 0.087960, Inf, 0.002247, 0.873605 +Loglikelihood: -61.702819 + +Status of colonist: 0, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -0.007835 +Status of colonist: 0, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: -0.008967 +Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 +Status of colonist: 4, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.497364 +Status of colonist: 2, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: 1.545883 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -13.576394 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.966530 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -8.752420 +Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.981045 +Parameters: +lambda^c, mu, K, gamma, lambda^a +0.195442, 0.087960, Inf, 0.002247, 0.873605 +Loglikelihood: -61.702819 + +Status of colonist: 0, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -0.007835 +Status of colonist: 0, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: -0.008967 +Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259847 +Status of colonist: 4, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.497364 +Status of colonist: 2, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: 1.545893 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -13.576392 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.966530 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -8.752420 +Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259847 +Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.981045 +Parameters: +lambda^c, mu, K, gamma, lambda^a +0.195442, 0.087960, Inf, 0.002247, 0.873605 +Loglikelihood: -61.702808 + +> +> +==2373== +==2373== HEAP SUMMARY: +==2373== in use at exit: 67,395,855 bytes in 14,084 blocks +==2373== total heap usage: 75,008,366 allocs, 74,994,282 frees, 40,220,411,885 bytes allocated +==2373== +==2373== LEAK SUMMARY: +==2373== definitely lost: 0 bytes in 0 blocks +==2373== indirectly lost: 0 bytes in 0 blocks +==2373== possibly lost: 0 bytes in 0 blocks +==2373== still reachable: 67,395,855 bytes in 14,084 blocks +==2373== of which reachable via heuristic: +==2373== newarray : 4,264 bytes in 1 blocks +==2373== suppressed: 0 bytes in 0 blocks +==2373== Reachable blocks (those to which a pointer was found) are not shown. +==2373== To see them, rerun with: --leak-check=full --show-leak-kinds=all +==2373== +==2373== For lists of detected and suppressed errors, rerun with: -s +==2373== ERROR SUMMARY: 0 errors from 0 contexts (suppressed: 0 from 0) diff --git a/src/DAISIE_loglik_rhs_FORTRAN.f95 b/src/DAISIE_loglik_rhs_FORTRAN.f95 index 2bb49994..5c121830 100644 --- a/src/DAISIE_loglik_rhs_FORTRAN.f95 +++ b/src/DAISIE_loglik_rhs_FORTRAN.f95 @@ -82,8 +82,8 @@ SUBROUTINE daisie_runmod (neq, t, Conc, dConc, yout, ip) INTEGER :: neq, ip(*), i, ii DOUBLE PRECISION :: t, Conc(2 * N + 1), dConc(2 * N + 1), yout(*) DOUBLE PRECISION :: xx1(N + 3), xx2(N + 3), xx3 - INTEGER :: il1(N), il2(N), il3in3(N), il4(N) - INTEGER :: in1(N), in2ix2(N) + INTEGER :: il1(N), il2(N), il3in3(N), il4(N) + INTEGER :: in1(N), in2ix2(N) INTEGER :: ix1(N), ix3(N), ix4(N) DOUBLE PRECISION :: laavec(N + 4 + 2 * kk),lacvec(N + 4 + 2 * kk) DOUBLE PRECISION :: muvec(N + 4 + 2 * kk),gamvec(N + 4 + 2 * kk) @@ -131,30 +131,30 @@ SUBROUTINE daisie_runmod (neq, t, Conc, dConc, yout, ip) xx1(1) = 0 xx1(2) = 0 xx2(1) = 0 - xx2(2) = 0 + xx2(2) = 0 DO I = 3, N + 2 - xx1(I) = Conc(I - 2) - xx2(I) = Conc(N + I - 2) - il1(I - 2) = I + kk - 1 - il2(I - 2) = I + kk + 1 - il3in3(I - 2) = I + kk - il4(I - 2) = I + kk - 2 - in1(I - 2) = I + 2 * kk - 1 - in2ix2(I - 2) = I + 1 - ix1(I - 2) = I - 1 - ix3(I - 2) = I - ix4(I - 2) = I - 2 - ENDDO - xx1(N + 3) = 0 - xx2(N + 3) = 0 + xx1(I) = Conc(I - 2) + xx2(I) = Conc(N + I - 2) + il1(I - 2) = I + kk - 1 + il2(I - 2) = I + kk + 1 + il3in3(I - 2) = I + kk + il4(I - 2) = I + kk - 2 + in1(I - 2) = I + 2 * kk - 1 + in2ix2(I - 2) = I + 1 + ix1(I - 2) = I - 1 + ix3(I - 2) = I + ix4(I - 2) = I - 2 + ENDDO + xx1(N + 3) = 0 + xx2(N + 3) = 0 xx3 = Conc(2 * N + 1) DO I = 1, N + 4 + 2 * kk - laavec(I) = P(I) - lacvec(I) = P(I + N + 4 + 2 * kk) - muvec(I) = P(I + 2 * (N + 4 + 2 * kk)) - gamvec(I) = P(I + 3 * (N + 4 + 2 * kk)) - nn(I) = P(I + 4 * (N + 4 + 2 * kk)) + laavec(I) = P(I) + lacvec(I) = P(I + N + 4 + 2 * kk) + muvec(I) = P(I + 2 * (N + 4 + 2 * kk)) + gamvec(I) = P(I + 3 * (N + 4 + 2 * kk)) + nn(I) = P(I + 4 * (N + 4 + 2 * kk)) ENDDO ! dx1 = laavec[il1 + 1] * xx2[ix1] + @@ -172,15 +172,15 @@ SUBROUTINE daisie_runmod (neq, t, Conc, dConc, yout, ip) ! dx3 <- 0 DO I = 1, N - dConc(I) = & - laavec(il1(I) + 1) * xx2(ix1(I)) + & - lacvec(il4(I) + 1) * xx2(ix4(I)) + & - muvec(il2(I) + 1) * xx2(ix3(I)) + & - lacvec(il1(I)) * nn(in1(I)) * xx1(ix1(I)) + & - muvec(il2(I)) * nn(in2ix2(I)) * xx1(in2ix2(I)) - & - (muvec(il3in3(I)) + lacvec(il3in3(I))) * & - nn(il3in3(I)) * xx1(ix3(I)) - & - gamvec(il3in3(I)) * xx1(ix3(I)) + dConc(I) = & + laavec(il1(I) + 1) * xx2(ix1(I)) + & + lacvec(il4(I) + 1) * xx2(ix4(I)) + & + muvec(il2(I) + 1) * xx2(ix3(I)) + & + lacvec(il1(I)) * nn(in1(I)) * xx1(ix1(I)) + & + muvec(il2(I)) * nn(in2ix2(I)) * xx1(in2ix2(I)) - & + (muvec(il3in3(I)) + lacvec(il3in3(I))) * & + nn(il3in3(I)) * xx1(ix3(I)) - & + gamvec(il3in3(I)) * xx1(ix3(I)) dConc(N + I) = & gamvec(il3in3(I)) * xx1(ix3(I)) + & lacvec(il1(I) + 1) * nn(in1(I)) * xx2(ix1(I)) + & @@ -188,7 +188,7 @@ SUBROUTINE daisie_runmod (neq, t, Conc, dConc, yout, ip) (muvec(il3in3(I) + 1) + lacvec(il3in3(I) + 1)) * & nn(il3in3(I) + 1) * xx2(ix3(I)) - & laavec(il3in3(I) + 1) * xx2(ix3(I)) - dConc(2 * N + 1) = 0 + dConc(2 * N + 1) = 0 ENDDO @@ -209,8 +209,8 @@ SUBROUTINE daisie_runmod1 (neq, t, Conc, dConc, yout, ip) INTEGER :: neq, ip(*), i, ii DOUBLE PRECISION :: t, Conc(4 * N), dConc(4 * N), yout(*) DOUBLE PRECISION :: xx1(N + 3), xx2(N + 3), xx3(N + 3), xx4(N + 3) - INTEGER :: il1(N), il2(N), il3in3(N), il4(N) - INTEGER :: in1(N), in2ix2(N), in4ix1(N) + INTEGER :: il1(N), il2(N), il3in3(N), il4(N) + INTEGER :: in1(N), in2ix2(N), in4ix1(N) INTEGER :: ix3(N), ix4(N) DOUBLE PRECISION :: laavec(N + 4 + 2 * kk),lacvec(N + 4 + 2 * kk) DOUBLE PRECISION :: muvec(N + 4 + 2 * kk),gamvec(N + 4 + 2 * kk) @@ -262,37 +262,37 @@ SUBROUTINE daisie_runmod1 (neq, t, Conc, dConc, yout, ip) xx1(1) = 0 xx1(2) = 0 xx2(1) = 0 - xx2(2) = 0 + xx2(2) = 0 xx3(1) = 0 - xx3(2) = 0 + xx3(2) = 0 xx4(1) = 0 - xx4(2) = 0 + xx4(2) = 0 DO I = 3, N + 2 - xx1(I) = Conc(I - 2) - xx2(I) = Conc(N + I - 2) - xx3(I) = Conc(2 * N + I - 2) - xx4(I) = Conc(3 * N + I - 2) - il1(I - 2) = I + kk - 1 - il2(I - 2) = I + kk + 1 - il3in3(I - 2) = I + kk - il4(I - 2) = I + kk - 2 - in1(I - 2) = I + 2 * kk - 1 - in2ix2(I - 2) = I + 1 - in4ix1(I - 2) = I - 1 - ix3(I - 2) = I - ix4(I - 2) = I - 2 - ENDDO - xx1(N + 3) = 0 - xx2(N + 3) = 0 + xx1(I) = Conc(I - 2) + xx2(I) = Conc(N + I - 2) + xx3(I) = Conc(2 * N + I - 2) + xx4(I) = Conc(3 * N + I - 2) + il1(I - 2) = I + kk - 1 + il2(I - 2) = I + kk + 1 + il3in3(I - 2) = I + kk + il4(I - 2) = I + kk - 2 + in1(I - 2) = I + 2 * kk - 1 + in2ix2(I - 2) = I + 1 + in4ix1(I - 2) = I - 1 + ix3(I - 2) = I + ix4(I - 2) = I - 2 + ENDDO + xx1(N + 3) = 0 + xx2(N + 3) = 0 xx3(N + 3) = 0 xx4(N + 3) = 0 DO I = 1, N + 4 + 2 * kk - laavec(I) = P(I) - lacvec(I) = P(I + N + 4 + 2 * kk) - muvec(I) = P(I + 2 * (N + 4 + 2 * kk)) - gamvec(I) = P(I + 3 * (N + 4 + 2 * kk)) - nn(I) = P(I + 4 * (N + 4 + 2 * kk)) + laavec(I) = P(I) + lacvec(I) = P(I + N + 4 + 2 * kk) + muvec(I) = P(I + 2 * (N + 4 + 2 * kk)) + gamvec(I) = P(I + 3 * (N + 4 + 2 * kk)) + nn(I) = P(I + 4 * (N + 4 + 2 * kk)) ENDDO ! dx1 <- lacvec[il1] * xx1[ix1] + @@ -304,15 +304,15 @@ SUBROUTINE daisie_runmod1 (neq, t, Conc, dConc, yout, ip) ! -gamvec[il3] * xx1[ix3] DO I = 1, N - dConc(I) = & - lacvec(il1(I)) * nn(in1(I)) * xx1(in4ix1(I)) + & - laavec(il1(I) + 1) * xx2(in4ix1(I)) + & - lacvec(il4(I) + 1) * xx2(ix4(I)) + & - muvec(il2(I)) * nn(in2ix2(I)) * xx1(in2ix2(I)) + & - muvec(il3in3(I) + 1) * xx2(ix3(I)) - & - (muvec(il3in3(I)) + lacvec(il3in3(I))) * & - nn(il3in3(I)) * xx1(ix3(I)) - & - gamvec(il3in3(I)) * xx1(ix3(I)) + dConc(I) = & + lacvec(il1(I)) * nn(in1(I)) * xx1(in4ix1(I)) + & + laavec(il1(I) + 1) * xx2(in4ix1(I)) + & + lacvec(il4(I) + 1) * xx2(ix4(I)) + & + muvec(il2(I)) * nn(in2ix2(I)) * xx1(in2ix2(I)) + & + muvec(il3in3(I) + 1) * xx2(ix3(I)) - & + (muvec(il3in3(I)) + lacvec(il3in3(I))) * & + nn(il3in3(I)) * xx1(ix3(I)) - & + gamvec(il3in3(I)) * xx1(ix3(I)) ! dx2 <- gamvec[il3] * xx1[ix3] + ! gamvec[il3] * xx3[ix3] + @@ -326,11 +326,11 @@ SUBROUTINE daisie_runmod1 (neq, t, Conc, dConc, yout, ip) gamvec(il3in3(I)) * xx1(ix3(I)) + & gamvec(il3in3(I)) * xx3(ix3(I)) + & gamvec(il3in3(I) + 1) * xx4(ix3(I)) + & - lacvec(il1(I) + 1) * nn(in1(I)) * xx2(in4ix1(I)) + & - muvec(il2(I) + 1) * nn(in2ix2(I)) * xx2(in2ix2(I)) - & - (muvec(il3in3(I) + 1) + lacvec(il3in3(I) + 1)) * & - nn(il3in3(I) + 1) * xx2(ix3(I)) - & - laavec(il3in3(I) + 1) * xx2(ix3(I)) + lacvec(il1(I) + 1) * nn(in1(I)) * xx2(in4ix1(I)) + & + muvec(il2(I) + 1) * nn(in2ix2(I)) * xx2(in2ix2(I)) - & + (muvec(il3in3(I) + 1) + lacvec(il3in3(I) + 1)) * & + nn(il3in3(I) + 1) * xx2(ix3(I)) - & + laavec(il3in3(I) + 1) * xx2(ix3(I)) ! dx3 <- lacvec[il1] * nn[in1] * xx3[ix1] + ! laavec[il1 + 1] * xx4[ix1] + @@ -355,7 +355,7 @@ SUBROUTINE daisie_runmod1 (neq, t, Conc, dConc, yout, ip) ! -(lacvec[il3 + 1] + muvec[il3 + 1]) * nn[in3 + 1] * xx4[ix3] + ! -gamvec[il3 + 1] * xx4[ix3] - dConc(3 * N + I) = & + dConc(3 * N + I) = & lacvec(il1(I) + 1) * nn(in1(I)) * xx4(in4ix1(I)) + & muvec(il2(I) + 1) * nn(in2ix2(I)) * xx4(in2ix2(I)) - & (lacvec(il3in3(I) + 1) + muvec(il3in3(I) + 1)) * & @@ -382,8 +382,8 @@ SUBROUTINE daisie_runmod2 (neq, t, Conc, dConc, yout, ip) INTEGER :: neq, ip(*), i, ii DOUBLE PRECISION :: t, Conc(3 * N), dConc(3 * N), yout(*) DOUBLE PRECISION :: xx1(N + 3), xx2(N + 3), xx3(N + 3) - INTEGER :: il1(N), il2(N), il3in3(N), il4(N) - INTEGER :: in1(N), in2ix2(N), in4ix1(N) + INTEGER :: il1(N), il2(N), il3in3(N), il4(N) + INTEGER :: in1(N), in2ix2(N), in4ix1(N) INTEGER :: ix3(N), ix4(N) DOUBLE PRECISION :: laavec(N + 4 + 2 * kk),lacvec(N + 4 + 2 * kk) DOUBLE PRECISION :: muvec(N + 4 + 2 * kk),gamvec(N + 4 + 2 * kk) @@ -433,33 +433,33 @@ SUBROUTINE daisie_runmod2 (neq, t, Conc, dConc, yout, ip) xx1(1) = 0 xx1(2) = 0 xx2(1) = 0 - xx2(2) = 0 + xx2(2) = 0 xx3(1) = 0 - xx3(2) = 0 + xx3(2) = 0 DO I = 3, N + 2 - xx1(I) = Conc(I - 2) - xx2(I) = Conc(N + I - 2) - xx3(I) = Conc(2 * N + I - 2) - il1(I - 2) = I + kk - 1 - il2(I - 2) = I + kk + 1 - il3in3(I - 2) = I + kk - il4(I - 2) = I + kk - 2 - in1(I - 2) = I + 2 * kk - 1 - in2ix2(I - 2) = I + 1 - in4ix1(I - 2) = I - 1 - ix3(I - 2) = I - ix4(I - 2) = I - 2 - ENDDO - xx1(N + 3) = 0 - xx2(N + 3) = 0 + xx1(I) = Conc(I - 2) + xx2(I) = Conc(N + I - 2) + xx3(I) = Conc(2 * N + I - 2) + il1(I - 2) = I + kk - 1 + il2(I - 2) = I + kk + 1 + il3in3(I - 2) = I + kk + il4(I - 2) = I + kk - 2 + in1(I - 2) = I + 2 * kk - 1 + in2ix2(I - 2) = I + 1 + in4ix1(I - 2) = I - 1 + ix3(I - 2) = I + ix4(I - 2) = I - 2 + ENDDO + xx1(N + 3) = 0 + xx2(N + 3) = 0 xx3(N + 3) = 0 DO I = 1, N + 4 + 2 * kk - laavec(I) = P(I) - lacvec(I) = P(I + N + 4 + 2 * kk) - muvec(I) = P(I + 2 * (N + 4 + 2 * kk)) - gamvec(I) = P(I + 3 * (N + 4 + 2 * kk)) - nn(I) = P(I + 4 * (N + 4 + 2 * kk)) + laavec(I) = P(I) + lacvec(I) = P(I + N + 4 + 2 * kk) + muvec(I) = P(I + 2 * (N + 4 + 2 * kk)) + gamvec(I) = P(I + 3 * (N + 4 + 2 * kk)) + nn(I) = P(I + 4 * (N + 4 + 2 * kk)) ENDDO ! dx1 = (laavec[il3] * xx3[ix3] + @@ -484,33 +484,33 @@ SUBROUTINE daisie_runmod2 (neq, t, Conc, dConc, yout, ip) ! -(laavec[il3] + gamvec[il3]) * xx3[ix3] DO I = 1, N - dConc(I) = 0 - IF(kk .EQ. 1) THEN - dConc(I) = laavec(il3in3(I)) * xx3(ix3(I)) + & - 2 * lacvec(il1(I)) * xx3(in4ix1(I)) - ENDIF - dConc(I) = dConc(I) + & - laavec(il1(I) + 1) * xx2(in4ix1(I)) + & - lacvec(il4(I) + 1) * xx2(ix4(I)) + & - muvec(il2(I) + 1) * xx2(ix3(I)) + & - lacvec(il1(I)) * nn(in1(I)) * xx1(in4ix1(I)) + & - muvec(il2(I)) * nn(in2ix2(I)) * xx1(in2ix2(I)) - & - (muvec(il3in3(I)) + lacvec(il3in3(I))) * & - nn(il3in3(I)) * xx1(ix3(I)) - & - gamvec(il3in3(I)) * xx1(ix3(I)) - dConc(N + I) = & - gamvec(il3in3(I)) * xx1(ix3(I)) + & - lacvec(il1(I) + 1) * nn(in1(I)) * xx2(in4ix1(I)) + & - muvec(il2(I) + 1) * nn(in2ix2(I)) * xx2(in2ix2(I)) - & - (muvec(il3in3(I) + 1) + lacvec(il3in3(I) + 1)) * & - nn(il3in3(I) + 1) * xx2(ix3(I)) - & - laavec(il3in3(I) + 1) * xx2(ix3(I)) - dConc(2 * N + I) = & - lacvec(il1(I)) * nn(in4ix1(I)) * xx3(in4ix1(I)) + & - muvec(il2(I)) * nn(in2ix2(I)) * xx3(in2ix2(I)) - & - (lacvec(il3in3(I)) + muvec(il3in3(I))) * & - nn(il3in3(I)) * xx3(ix3(I)) - & - (laavec(il3in3(I)) + gamvec(il3in3(I))) * xx3(ix3(I)) + dConc(I) = 0 + IF(kk .EQ. 1) THEN + dConc(I) = laavec(il3in3(I)) * xx3(ix3(I)) + & + 2 * lacvec(il1(I)) * xx3(in4ix1(I)) + ENDIF + dConc(I) = dConc(I) + & + laavec(il1(I) + 1) * xx2(in4ix1(I)) + & + lacvec(il4(I) + 1) * xx2(ix4(I)) + & + muvec(il2(I) + 1) * xx2(ix3(I)) + & + lacvec(il1(I)) * nn(in1(I)) * xx1(in4ix1(I)) + & + muvec(il2(I)) * nn(in2ix2(I)) * xx1(in2ix2(I)) - & + (muvec(il3in3(I)) + lacvec(il3in3(I))) * & + nn(il3in3(I)) * xx1(ix3(I)) - & + gamvec(il3in3(I)) * xx1(ix3(I)) + dConc(N + I) = & + gamvec(il3in3(I)) * xx1(ix3(I)) + & + lacvec(il1(I) + 1) * nn(in1(I)) * xx2(in4ix1(I)) + & + muvec(il2(I) + 1) * nn(in2ix2(I)) * xx2(in2ix2(I)) - & + (muvec(il3in3(I) + 1) + lacvec(il3in3(I) + 1)) * & + nn(il3in3(I) + 1) * xx2(ix3(I)) - & + laavec(il3in3(I) + 1) * xx2(ix3(I)) + dConc(2 * N + I) = & + lacvec(il1(I)) * nn(in4ix1(I)) * xx3(in4ix1(I)) + & + muvec(il2(I)) * nn(in2ix2(I)) * xx3(in2ix2(I)) - & + (lacvec(il3in3(I)) + muvec(il3in3(I))) * & + nn(il3in3(I)) * xx3(ix3(I)) - & + (laavec(il3in3(I)) + gamvec(il3in3(I))) * xx3(ix3(I)) ENDDO END SUBROUTINE daisie_runmod2 diff --git a/src/DAISIE_odeint.h b/src/DAISIE_odeint.h index bb84dd62..ec63abde 100644 --- a/src/DAISIE_odeint.h +++ b/src/DAISIE_odeint.h @@ -6,6 +6,7 @@ #include "config.h" #include "ublas_types.h" +#include "patched_bulrisch_stoer.hpp" // shadow buggy boost header #include #include #include diff --git a/src/config.h b/src/config.h index a14c60c8..32cfa85b 100644 --- a/src/config.h +++ b/src/config.h @@ -8,18 +8,21 @@ // moved to Makevars[.win] /* #if __cplusplus >= 201703L -#ifdef _HAS_AUTO_PTR_ETC -#undef _HAS_AUTO_PTR_ETC -#endif -#define _HAS_AUTO_PTR_ETC 0 +# ifdef _HAS_AUTO_PTR_ETC +# undef _HAS_AUTO_PTR_ETC +# endif +# define _HAS_AUTO_PTR_ETC 0 #endif */ // Special case to make use of some steppers that would include // boost/get_pointer.hpp -#ifdef BOOST_NO_AUTO_PTR -#undef BOOST_NO_AUTO_PTR +#ifndef BOOST_NO_AUTO_PTR +# define BOOST_NO_AUTO_PTR #endif -#define BOOST_NO_AUTO_PTR + +// uncomment if unitialized member variable bulirsch_stoer::m_dt_last +// is fixed in boost (BH) +#define USE_BULRISCH_STOER_PATCH #endif diff --git a/src/patched_bulrisch_stoer.hpp b/src/patched_bulrisch_stoer.hpp new file mode 100644 index 00000000..cad4f521 --- /dev/null +++ b/src/patched_bulrisch_stoer.hpp @@ -0,0 +1,661 @@ +/* + Patched version of + boost/numeric/odeint/stepper/bulirsch_stoer.hpp + + Addresses unitialized member variable bulirsch_stoer<>::m_dt_last. + + Include this header before boost/numeric/odeint to shadow + boost/numeriuc/odeint/stepper/bulrisch_stoer.hpp + + The issue is *not* fixed in BOOST_VERSION == 1.81.0. + We need to check for fixes in upcomming boost (BH) releases. + + Hanno Hildenbrandt 2023 +*/ + +/* + [auto_generated] + boost/numeric/odeint/stepper/bulirsch_stoer.hpp + + [begin_description] + Implementation of the Burlish-Stoer method. As described in + Ernst Hairer, Syvert Paul Norsett, Gerhard Wanner + Solving Ordinary Differential Equations I. Nonstiff Problems. + Springer Series in Comput. Mathematics, Vol. 8, Springer-Verlag 1987, Second revised edition 1993. + [end_description] + + Copyright 2011-2013 Mario Mulansky + Copyright 2011-2013 Karsten Ahnert + Copyright 2012 Christoph Koke + + Distributed under the Boost Software License, Version 1.0. + (See accompanying file LICENSE_1_0.txt or + copy at http://www.boost.org/LICENSE_1_0.txt) +*/ + +#ifdef USE_BULRISCH_STOER_PATCH + +#ifndef BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_HPP_INCLUDED +#define BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_HPP_INCLUDED + + +#include + +#include + +#include // for min/max guidelines + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace boost { +namespace numeric { +namespace odeint { + +template< + class State , + class Value = double , + class Deriv = State , + class Time = Value , + class Algebra = typename algebra_dispatcher< State >::algebra_type , + class Operations = typename operations_dispatcher< State >::operations_type , + class Resizer = initially_resizer + > +class bulirsch_stoer { + +public: + + typedef State state_type; + typedef Value value_type; + typedef Deriv deriv_type; + typedef Time time_type; + typedef Algebra algebra_type; + typedef Operations operations_type; + typedef Resizer resizer_type; +#ifndef DOXYGEN_SKIP + typedef state_wrapper< state_type > wrapped_state_type; + typedef state_wrapper< deriv_type > wrapped_deriv_type; + typedef controlled_stepper_tag stepper_category; + + typedef bulirsch_stoer< State , Value , Deriv , Time , Algebra , Operations , Resizer > controlled_error_bs_type; + + typedef typename inverse_time< time_type >::type inv_time_type; + + typedef std::vector< value_type > value_vector; + typedef std::vector< time_type > time_vector; + typedef std::vector< inv_time_type > inv_time_vector; //should be 1/time_type for boost.units + typedef std::vector< value_vector > value_matrix; + typedef std::vector< size_t > int_vector; + typedef std::vector< wrapped_state_type > state_table_type; +#endif //DOXYGEN_SKIP + const static size_t m_k_max = 8; + + bulirsch_stoer( + value_type eps_abs = 1E-6 , value_type eps_rel = 1E-6 , + value_type factor_x = 1.0 , value_type factor_dxdt = 1.0 , + time_type max_dt = static_cast(0)) + : m_error_checker( eps_abs , eps_rel , factor_x, factor_dxdt ) , m_midpoint() , + m_last_step_rejected( false ) , m_first( true ) , + m_dt_last(static_cast(-1)) , // !!! patched !!! + m_max_dt(max_dt) , + m_interval_sequence( m_k_max+1 ) , + m_coeff( m_k_max+1 ) , + m_cost( m_k_max+1 ) , + m_facmin_table( m_k_max+1 ) , + m_table( m_k_max ) , + STEPFAC1( 0.65 ) , STEPFAC2( 0.94 ) , STEPFAC3( 0.02 ) , STEPFAC4( 4.0 ) , KFAC1( 0.8 ) , KFAC2( 0.9 ) + { + BOOST_USING_STD_MIN(); + BOOST_USING_STD_MAX(); + /* initialize sequence of stage numbers and work */ + for( unsigned short i = 0; i < m_k_max+1; i++ ) + { + m_interval_sequence[i] = 2 * (i+1); + if( i == 0 ) + m_cost[i] = m_interval_sequence[i]; + else + m_cost[i] = m_cost[i-1] + m_interval_sequence[i]; + m_coeff[i].resize(i); + m_facmin_table[i] = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , static_cast< value_type >(1) / static_cast< value_type >( 2*i+1 ) ); + for( size_t k = 0 ; k < i ; ++k ) + { + const value_type r = static_cast< value_type >( m_interval_sequence[i] ) / static_cast< value_type >( m_interval_sequence[k] ); + m_coeff[i][k] = 1.0 / ( r*r - static_cast< value_type >( 1.0 ) ); // coefficients for extrapolation + } + } + reset(); + } + + + /* + * Version 1 : try_step( sys , x , t , dt ) + * + * The overloads are needed to solve the forwarding problem + */ + template< class System , class StateInOut > + controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt ) + { + return try_step_v1( system , x , t, dt ); + } + + /** + * \brief Second version to solve the forwarding problem, can be used with Boost.Range as StateInOut. + */ + template< class System , class StateInOut > + controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt ) + { + return try_step_v1( system , x , t, dt ); + } + + /* + * Version 2 : try_step( sys , x , dxdt , t , dt ) + * + * this version does not solve the forwarding problem, boost.range can not be used + */ + template< class System , class StateInOut , class DerivIn > + controlled_step_result try_step( System system , StateInOut &x , const DerivIn &dxdt , time_type &t , time_type &dt ) + { + m_xnew_resizer.adjust_size( x , detail::bind( &controlled_error_bs_type::template resize_m_xnew< StateInOut > , detail::ref( *this ) , detail::_1 ) ); + controlled_step_result res = try_step( system , x , dxdt , t , m_xnew.m_v , dt ); + if( res == success ) + { + boost::numeric::odeint::copy( m_xnew.m_v , x ); + } + return res; + } + + /* + * Version 3 : try_step( sys , in , t , out , dt ) + * + * this version does not solve the forwarding problem, boost.range can not be used + */ + template< class System , class StateIn , class StateOut > + typename boost::disable_if< boost::is_same< StateIn , time_type > , controlled_step_result >::type + try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt ) + { + typename odeint::unwrap_reference< System >::type &sys = system; + m_dxdt_resizer.adjust_size( in , detail::bind( &controlled_error_bs_type::template resize_m_dxdt< StateIn > , detail::ref( *this ) , detail::_1 ) ); + sys( in , m_dxdt.m_v , t ); + return try_step( system , in , m_dxdt.m_v , t , out , dt ); + } + + + /* + * Full version : try_step( sys , in , dxdt_in , t , out , dt ) + * + * contains the actual implementation + */ + template< class System , class StateIn , class DerivIn , class StateOut > + controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt ) + { + if( m_max_dt != static_cast(0) && detail::less_with_sign(m_max_dt, dt, dt) ) + { + // given step size is bigger then max_dt + // set limit and return fail + dt = m_max_dt; + return fail; + } + + BOOST_USING_STD_MIN(); + BOOST_USING_STD_MAX(); + + static const value_type val1( 1.0 ); + + if( m_resizer.adjust_size( in , detail::bind( &controlled_error_bs_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ) ) + { + reset(); // system resized -> reset + } + + if( dt != m_dt_last ) + { + reset(); // step size changed from outside -> reset + } + + bool reject( true ); + + time_vector h_opt( m_k_max+1 ); + inv_time_vector work( m_k_max+1 ); + + time_type new_h = dt; + + /* m_current_k_opt is the estimated current optimal stage number */ + for( size_t k = 0 ; k <= m_current_k_opt+1 ; k++ ) + { + /* the stage counts are stored in m_interval_sequence */ + m_midpoint.set_steps( m_interval_sequence[k] ); + if( k == 0 ) + { + m_midpoint.do_step( system , in , dxdt , t , out , dt ); + /* the first step, nothing more to do */ + } + else + { + m_midpoint.do_step( system , in , dxdt , t , m_table[k-1].m_v , dt ); + extrapolate( k , m_table , m_coeff , out ); + // get error estimate + m_algebra.for_each3( m_err.m_v , out , m_table[0].m_v , + typename operations_type::template scale_sum2< value_type , value_type >( val1 , -val1 ) ); + const value_type error = m_error_checker.error( m_algebra , in , dxdt , m_err.m_v , dt ); + h_opt[k] = calc_h_opt( dt , error , k ); + work[k] = static_cast( m_cost[k] ) / h_opt[k]; + + if( (k == m_current_k_opt-1) || m_first ) + { // convergence before k_opt ? + if( error < 1.0 ) + { + //convergence + reject = false; + if( (work[k] < KFAC2*work[k-1]) || (m_current_k_opt <= 2) ) + { + // leave order as is (except we were in first round) + m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast(m_k_max)-1 , max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast(k)+1 ) ); + new_h = h_opt[k]; + new_h *= static_cast( m_cost[k+1] ) / static_cast( m_cost[k] ); + } else { + m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast(m_k_max)-1 , max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast(k) ) ); + new_h = h_opt[k]; + } + break; + } + else if( should_reject( error , k ) && !m_first ) + { + reject = true; + new_h = h_opt[k]; + break; + } + } + if( k == m_current_k_opt ) + { // convergence at k_opt ? + if( error < 1.0 ) + { + //convergence + reject = false; + if( (work[k-1] < KFAC2*work[k]) ) + { + m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast(m_current_k_opt)-1 ); + new_h = h_opt[m_current_k_opt]; + } + else if( (work[k] < KFAC2*work[k-1]) && !m_last_step_rejected ) + { + m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast(m_k_max-1) , static_cast(m_current_k_opt)+1 ); + new_h = h_opt[k]; + new_h *= static_cast(m_cost[m_current_k_opt])/static_cast(m_cost[k]); + } else + new_h = h_opt[m_current_k_opt]; + break; + } + else if( should_reject( error , k ) ) + { + reject = true; + new_h = h_opt[m_current_k_opt]; + break; + } + } + if( k == m_current_k_opt+1 ) + { // convergence at k_opt+1 ? + if( error < 1.0 ) + { //convergence + reject = false; + if( work[k-2] < KFAC2*work[k-1] ) + m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast(m_current_k_opt)-1 ); + if( (work[k] < KFAC2*work[m_current_k_opt]) && !m_last_step_rejected ) + m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast(m_k_max)-1 , static_cast(k) ); + new_h = h_opt[m_current_k_opt]; + } else + { + reject = true; + new_h = h_opt[m_current_k_opt]; + } + break; + } + } + } + + if( !reject ) + { + t += dt; + } + + if( !m_last_step_rejected || boost::numeric::odeint::detail::less_with_sign(new_h, dt, dt) ) + { + // limit step size + if( m_max_dt != static_cast(0) ) + { + new_h = detail::min_abs(m_max_dt, new_h); + } + m_dt_last = new_h; + dt = new_h; + } + + m_last_step_rejected = reject; + m_first = false; + + if( reject ) + return fail; + else + return success; + } + + /** \brief Resets the internal state of the stepper */ + void reset() + { + m_first = true; + m_last_step_rejected = false; + // crude estimate of optimal order + m_current_k_opt = 4; + /* no calculation because log10 might not exist for value_type! + const value_type logfact( -log10( max BOOST_PREVENT_MACRO_SUBSTITUTION( eps_rel , static_cast< value_type >(1.0E-12) ) ) * 0.6 + 0.5 ); + m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast( 1 ) , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast( m_k_max-1 ) , logfact )); + */ + } + + + /* Resizer methods */ + + template< class StateIn > + void adjust_size( const StateIn &x ) + { + resize_m_dxdt( x ); + resize_m_xnew( x ); + resize_impl( x ); + m_midpoint.adjust_size( x ); + } + + +private: + + template< class StateIn > + bool resize_m_dxdt( const StateIn &x ) + { + return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable::type() ); + } + + template< class StateIn > + bool resize_m_xnew( const StateIn &x ) + { + return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable::type() ); + } + + template< class StateIn > + bool resize_impl( const StateIn &x ) + { + bool resized( false ); + for( size_t i = 0 ; i < m_k_max ; ++i ) + resized |= adjust_size_by_resizeability( m_table[i] , x , typename is_resizeable::type() ); + resized |= adjust_size_by_resizeability( m_err , x , typename is_resizeable::type() ); + return resized; + } + + + template< class System , class StateInOut > + controlled_step_result try_step_v1( System system , StateInOut &x , time_type &t , time_type &dt ) + { + typename odeint::unwrap_reference< System >::type &sys = system; + m_dxdt_resizer.adjust_size( x , detail::bind( &controlled_error_bs_type::template resize_m_dxdt< StateInOut > , detail::ref( *this ) , detail::_1 ) ); + sys( x , m_dxdt.m_v ,t ); + return try_step( system , x , m_dxdt.m_v , t , dt ); + } + + + template< class StateInOut > + void extrapolate( size_t k , state_table_type &table , const value_matrix &coeff , StateInOut &xest ) + /* polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf + uses the obtained intermediate results to extrapolate to dt->0 + */ + { + static const value_type val1 = static_cast< value_type >( 1.0 ); + for( int j=k-1 ; j>0 ; --j ) + { + m_algebra.for_each3( table[j-1].m_v , table[j].m_v , table[j-1].m_v , + typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k][j] , -coeff[k][j] ) ); + } + m_algebra.for_each3( xest , table[0].m_v , xest , + typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k][0] , -coeff[k][0]) ); + } + + time_type calc_h_opt( time_type h , value_type error , size_t k ) const + /* calculates the optimal step size for a given error and stage number */ + { + BOOST_USING_STD_MIN(); + BOOST_USING_STD_MAX(); + using std::pow; + value_type expo( 1.0/(2*k+1) ); + value_type facmin = m_facmin_table[k]; + value_type fac; + if (error == 0.0) + fac=1.0/facmin; + else + { + fac = STEPFAC2 / pow BOOST_PREVENT_MACRO_SUBSTITUTION( error / STEPFAC1 , expo ); + fac = max BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast(facmin/STEPFAC4) , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast(1.0/facmin) , fac ) ); + } + return h*fac; + } + + controlled_step_result set_k_opt( size_t k , const inv_time_vector &work , const time_vector &h_opt , time_type &dt ) + /* calculates the optimal stage number */ + { + if( k == 1 ) + { + m_current_k_opt = 2; + return success; + } + if( (work[k-1] < KFAC1*work[k]) || (k == m_k_max) ) + { // order decrease + m_current_k_opt = k-1; + dt = h_opt[ m_current_k_opt ]; + return success; + } + else if( (work[k] < KFAC2*work[k-1]) || m_last_step_rejected || (k == m_k_max-1) ) + { // same order - also do this if last step got rejected + m_current_k_opt = k; + dt = h_opt[ m_current_k_opt ]; + return success; + } + else + { // order increase - only if last step was not rejected + m_current_k_opt = k+1; + dt = h_opt[ m_current_k_opt-1 ] * m_cost[ m_current_k_opt ] / m_cost[ m_current_k_opt-1 ] ; + return success; + } + } + + bool in_convergence_window( size_t k ) const + { + if( (k == m_current_k_opt-1) && !m_last_step_rejected ) + return true; // decrease stepsize only if last step was not rejected + return ( (k == m_current_k_opt) || (k == m_current_k_opt+1) ); + } + + bool should_reject( value_type error , size_t k ) const + { + if( k == m_current_k_opt-1 ) + { + const value_type d = m_interval_sequence[m_current_k_opt] * m_interval_sequence[m_current_k_opt+1] / + (m_interval_sequence[0]*m_interval_sequence[0]); + //step will fail, criterion 17.3.17 in NR + return ( error > d*d ); + } + else if( k == m_current_k_opt ) + { + const value_type d = m_interval_sequence[m_current_k_opt] / m_interval_sequence[0]; + return ( error > d*d ); + } else + return error > 1.0; + } + + default_error_checker< value_type, algebra_type , operations_type > m_error_checker; + modified_midpoint< state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > m_midpoint; + + bool m_last_step_rejected; + bool m_first; + + time_type m_dt_last; + time_type m_t_last; + time_type m_max_dt; + + size_t m_current_k_opt; + + algebra_type m_algebra; + + resizer_type m_dxdt_resizer; + resizer_type m_xnew_resizer; + resizer_type m_resizer; + + wrapped_state_type m_xnew; + wrapped_state_type m_err; + wrapped_deriv_type m_dxdt; + + int_vector m_interval_sequence; // stores the successive interval counts + value_matrix m_coeff; + int_vector m_cost; // costs for interval count + value_vector m_facmin_table; // for precomputed facmin to save pow calls + + state_table_type m_table; // sequence of states for extrapolation + + value_type STEPFAC1 , STEPFAC2 , STEPFAC3 , STEPFAC4 , KFAC1 , KFAC2; +}; + + +/******** DOXYGEN ********/ +/** + * \class bulirsch_stoer + * \brief The Bulirsch-Stoer algorithm. + * + * The Bulirsch-Stoer is a controlled stepper that adjusts both step size + * and order of the method. The algorithm uses the modified midpoint and + * a polynomial extrapolation compute the solution. + * + * \tparam State The state type. + * \tparam Value The value type. + * \tparam Deriv The type representing the time derivative of the state. + * \tparam Time The time representing the independent variable - the time. + * \tparam Algebra The algebra type. + * \tparam Operations The operations type. + * \tparam Resizer The resizer policy type. + */ + + /** + * \fn bulirsch_stoer::bulirsch_stoer( value_type eps_abs , value_type eps_rel , value_type factor_x , value_type factor_dxdt ) + * \brief Constructs the bulirsch_stoer class, including initialization of + * the error bounds. + * + * \param eps_abs Absolute tolerance level. + * \param eps_rel Relative tolerance level. + * \param factor_x Factor for the weight of the state. + * \param factor_dxdt Factor for the weight of the derivative. + */ + + /** + * \fn bulirsch_stoer::try_step( System system , StateInOut &x , time_type &t , time_type &dt ) + * \brief Tries to perform one step. + * + * This method tries to do one step with step size dt. If the error estimate + * is to large, the step is rejected and the method returns fail and the + * step size dt is reduced. If the error estimate is acceptably small, the + * step is performed, success is returned and dt might be increased to make + * the steps as large as possible. This method also updates t if a step is + * performed. Also, the internal order of the stepper is adjusted if required. + * + * \param system The system function to solve, hence the r.h.s. of the ODE. + * It must fulfill the Simple System concept. + * \param x The state of the ODE which should be solved. Overwritten if + * the step is successful. + * \param t The value of the time. Updated if the step is successful. + * \param dt The step size. Updated. + * \return success if the step was accepted, fail otherwise. + */ + + /** + * \fn bulirsch_stoer::try_step( System system , StateInOut &x , const DerivIn &dxdt , time_type &t , time_type &dt ) + * \brief Tries to perform one step. + * + * This method tries to do one step with step size dt. If the error estimate + * is to large, the step is rejected and the method returns fail and the + * step size dt is reduced. If the error estimate is acceptably small, the + * step is performed, success is returned and dt might be increased to make + * the steps as large as possible. This method also updates t if a step is + * performed. Also, the internal order of the stepper is adjusted if required. + * + * \param system The system function to solve, hence the r.h.s. of the ODE. + * It must fulfill the Simple System concept. + * \param x The state of the ODE which should be solved. Overwritten if + * the step is successful. + * \param dxdt The derivative of state. + * \param t The value of the time. Updated if the step is successful. + * \param dt The step size. Updated. + * \return success if the step was accepted, fail otherwise. + */ + + /** + * \fn bulirsch_stoer::try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt ) + * \brief Tries to perform one step. + * + * \note This method is disabled if state_type=time_type to avoid ambiguity. + * + * This method tries to do one step with step size dt. If the error estimate + * is to large, the step is rejected and the method returns fail and the + * step size dt is reduced. If the error estimate is acceptably small, the + * step is performed, success is returned and dt might be increased to make + * the steps as large as possible. This method also updates t if a step is + * performed. Also, the internal order of the stepper is adjusted if required. + * + * \param system The system function to solve, hence the r.h.s. of the ODE. + * It must fulfill the Simple System concept. + * \param in The state of the ODE which should be solved. + * \param t The value of the time. Updated if the step is successful. + * \param out Used to store the result of the step. + * \param dt The step size. Updated. + * \return success if the step was accepted, fail otherwise. + */ + + + /** + * \fn bulirsch_stoer::try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt ) + * \brief Tries to perform one step. + * + * This method tries to do one step with step size dt. If the error estimate + * is to large, the step is rejected and the method returns fail and the + * step size dt is reduced. If the error estimate is acceptably small, the + * step is performed, success is returned and dt might be increased to make + * the steps as large as possible. This method also updates t if a step is + * performed. Also, the internal order of the stepper is adjusted if required. + * + * \param system The system function to solve, hence the r.h.s. of the ODE. + * It must fulfill the Simple System concept. + * \param in The state of the ODE which should be solved. + * \param dxdt The derivative of state. + * \param t The value of the time. Updated if the step is successful. + * \param out Used to store the result of the step. + * \param dt The step size. Updated. + * \return success if the step was accepted, fail otherwise. + */ + + + /** + * \fn bulirsch_stoer::adjust_size( const StateIn &x ) + * \brief Adjust the size of all temporaries in the stepper manually. + * \param x A state from which the size of the temporaries to be resized is deduced. + */ + +} +} +} + +#endif // BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_HPP_INCLUDED + +#endif // USE_BULRISCH_STOER_PATCH From d8a1293f94f33b67564845687c7a0a169aaf8e96 Mon Sep 17 00:00:00 2001 From: Hanno Hildenbrandt Date: Thu, 9 Mar 2023 12:58:45 +0100 Subject: [PATCH 04/13] nop --- source_odeint.R | 195 ------------------------------------------------ 1 file changed, 195 deletions(-) diff --git a/source_odeint.R b/source_odeint.R index 8a6ea202..5ec7ac01 100644 --- a/source_odeint.R +++ b/source_odeint.R @@ -45,198 +45,3 @@ methode <- 'odeint::adams_bashforth_moulton_1' DAISIE_CS_max_steps(100000000) DAISIE_abm_factor(0.000001) loglik_abm <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) - - -==2373== Memcheck, a memory error detector -==2373== Copyright (C) 2002-2022, and GNU GPL'd, by Julian Seward et al. -==2373== Using Valgrind-3.19.0 and LibVEX; rerun with -h for copyright info -==2373== Command: /home/hanno/opt/Rval/lib/R/bin/exec/R --vanilla -e source('source_odeint.R') -==2373== - -R Under development (unstable) (2023-03-07 r83950) -- "Unsuffered Consequences" -Copyright (C) 2023 The R Foundation for Statistical Computing -Platform: x86_64-pc-linux-gnu (64-bit) - -R is free software and comes with ABSOLUTELY NO WARRANTY. -You are welcome to redistribute it under certain conditions. -Type 'license()' or 'licence()' for distribution details. - -R is a collaborative project with many contributors. -Type 'contributors()' for more information and -'citation()' on how to cite R or R packages in publications. - -Type 'demo()' for some demos, 'help()' for on-line help, or -'help.start()' for an HTML browser interface to help. -Type 'q()' to quit R. - -> source('source_odeint.R') -Status of colonist: 0, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -0.007835 -Status of colonist: 0, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: -0.008967 -Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 -Status of colonist: 4, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.497364 -Status of colonist: 2, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: 1.545883 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -13.576394 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.966530 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -8.752420 -Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.981045 -Parameters: -lambda^c, mu, K, gamma, lambda^a -0.195442, 0.087960, Inf, 0.002247, 0.873605 -Loglikelihood: -61.702819 - -Status of colonist: 0, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -0.007835 -Status of colonist: 0, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: -0.008967 -Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 -Status of colonist: 4, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.497364 -Status of colonist: 2, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: 1.545883 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -13.576394 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.966530 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -8.752420 -Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.981045 -Parameters: -lambda^c, mu, K, gamma, lambda^a -0.195442, 0.087960, Inf, 0.002247, 0.873605 -Loglikelihood: -61.702819 - -Status of colonist: 0, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -0.007835 -Status of colonist: 0, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: -0.008967 -Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 -Status of colonist: 4, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.497364 -Status of colonist: 2, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: 1.545883 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -13.576394 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.966530 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -8.752420 -Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.981045 -Parameters: -lambda^c, mu, K, gamma, lambda^a -0.195442, 0.087960, Inf, 0.002247, 0.873605 -Loglikelihood: -61.702819 - -Status of colonist: 0, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -0.007835 -Status of colonist: 0, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: -0.008967 -Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 -Status of colonist: 4, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.497364 -Status of colonist: 2, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: 1.545883 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -13.576394 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.966530 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -8.752420 -Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.981045 -Parameters: -lambda^c, mu, K, gamma, lambda^a -0.195442, 0.087960, Inf, 0.002247, 0.873605 -Loglikelihood: -61.702819 - -Status of colonist: 0, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -0.007835 -Status of colonist: 0, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: -0.008967 -Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 -Status of colonist: 4, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.497364 -Status of colonist: 2, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: 1.545883 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -13.576394 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.966530 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -8.752420 -Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.981045 -Parameters: -lambda^c, mu, K, gamma, lambda^a -0.195442, 0.087960, Inf, 0.002247, 0.873605 -Loglikelihood: -61.702819 - -Status of colonist: 0, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -0.007835 -Status of colonist: 0, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: -0.008967 -Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 -Status of colonist: 4, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.497364 -Status of colonist: 2, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: 1.545883 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -13.576394 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.966530 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -8.752420 -Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.981045 -Parameters: -lambda^c, mu, K, gamma, lambda^a -0.195442, 0.087960, Inf, 0.002247, 0.873605 -Loglikelihood: -61.702819 - -Status of colonist: 0, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -0.007835 -Status of colonist: 0, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: -0.008967 -Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 -Status of colonist: 4, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.497364 -Status of colonist: 2, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: 1.545883 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -13.576394 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.966530 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -8.752420 -Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.981045 -Parameters: -lambda^c, mu, K, gamma, lambda^a -0.195442, 0.087960, Inf, 0.002247, 0.873605 -Loglikelihood: -61.702819 - -Status of colonist: 0, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -0.007835 -Status of colonist: 0, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: -0.008967 -Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 -Status of colonist: 4, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.497364 -Status of colonist: 2, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: 1.545883 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -13.576394 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.966530 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -8.752420 -Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.981045 -Parameters: -lambda^c, mu, K, gamma, lambda^a -0.195442, 0.087960, Inf, 0.002247, 0.873605 -Loglikelihood: -61.702819 - -Status of colonist: 0, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -0.007835 -Status of colonist: 0, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: -0.008967 -Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 -Status of colonist: 4, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.497364 -Status of colonist: 2, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: 1.545883 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -13.576394 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.966530 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -8.752420 -Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259846 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.981045 -Parameters: -lambda^c, mu, K, gamma, lambda^a -0.195442, 0.087960, Inf, 0.002247, 0.873605 -Loglikelihood: -61.702819 - -Status of colonist: 0, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -0.007835 -Status of colonist: 0, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: -0.008967 -Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259847 -Status of colonist: 4, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.497364 -Status of colonist: 2, Parameters: 3755.202241, 8.909285, 14.999999, 0.002247, 0.873605, Loglikelihood: 1.545893 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -13.576392 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.966530 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -8.752420 -Status of colonist: 1, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.259847 -Status of colonist: 2, Parameters: 0.195442, 0.087960, Inf, 0.002247, 0.873605, Loglikelihood: -6.981045 -Parameters: -lambda^c, mu, K, gamma, lambda^a -0.195442, 0.087960, Inf, 0.002247, 0.873605 -Loglikelihood: -61.702808 - -> -> -==2373== -==2373== HEAP SUMMARY: -==2373== in use at exit: 67,395,855 bytes in 14,084 blocks -==2373== total heap usage: 75,008,366 allocs, 74,994,282 frees, 40,220,411,885 bytes allocated -==2373== -==2373== LEAK SUMMARY: -==2373== definitely lost: 0 bytes in 0 blocks -==2373== indirectly lost: 0 bytes in 0 blocks -==2373== possibly lost: 0 bytes in 0 blocks -==2373== still reachable: 67,395,855 bytes in 14,084 blocks -==2373== of which reachable via heuristic: -==2373== newarray : 4,264 bytes in 1 blocks -==2373== suppressed: 0 bytes in 0 blocks -==2373== Reachable blocks (those to which a pointer was found) are not shown. -==2373== To see them, rerun with: --leak-check=full --show-leak-kinds=all -==2373== -==2373== For lists of detected and suppressed errors, rerun with: -s -==2373== ERROR SUMMARY: 0 errors from 0 contexts (suppressed: 0 from 0) From 52b008dd848f5ed3178dba844b673613397a2721 Mon Sep 17 00:00:00 2001 From: Neves-P Date: Thu, 9 Mar 2023 16:38:36 +0100 Subject: [PATCH 05/13] Remove TB and ABM1 stepper for separate test --- tests/testthat/{test_odeint.R => test-odeint.R} | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) rename tests/testthat/{test_odeint.R => test-odeint.R} (85%) diff --git a/tests/testthat/test_odeint.R b/tests/testthat/test-odeint.R similarity index 85% rename from tests/testthat/test_odeint.R rename to tests/testthat/test-odeint.R index 6a47f4a9..eb03992a 100644 --- a/tests/testthat/test_odeint.R +++ b/tests/testthat/test-odeint.R @@ -33,18 +33,11 @@ test_that("odeint solvers give the same result as deSolve solvers", { loglik_rkd5 <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) methode <- 'odeint::bulirsch_stoer' loglik_bs <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) - methode <- 'odeint::rosenbrock4' - loglik_rb <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) - methode <- 'odeint::adams_bashforth_moulton_1' - DAISIE_CS_max_steps(100000000) - DAISIE_abm_factor(0.000001) - loglik_abm <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) expect_equal(loglik_lsodes, loglik_rkck54) expect_equal(loglik_lsodes, loglik_rkf78) expect_equal(loglik_lsodes, loglik_rkd5) expect_equal(loglik_lsodes, loglik_bs) - expect_equal(loglik_lsodes, loglik_rb) - expect_equal(loglik_lsodes, loglik_abm, tol = 1E-6) + pars1a <- pars1 pars1a[6] <- Inf From ffee1ae7e28bbab73c6da84b13e12000b11a8858 Mon Sep 17 00:00:00 2001 From: Neves-P Date: Thu, 9 Mar 2023 16:39:26 +0100 Subject: [PATCH 06/13] Add (skipped) RB and ABM1 stepper tests These steppers take too long to test properly with automation. These tests are kept so they can be run manually. --- .../testthat/test-rosenbrock_abm1_steppers.R | 28 +++++++++++++++++++ 1 file changed, 28 insertions(+) create mode 100644 tests/testthat/test-rosenbrock_abm1_steppers.R diff --git a/tests/testthat/test-rosenbrock_abm1_steppers.R b/tests/testthat/test-rosenbrock_abm1_steppers.R new file mode 100644 index 00000000..6037871c --- /dev/null +++ b/tests/testthat/test-rosenbrock_abm1_steppers.R @@ -0,0 +1,28 @@ +test_that("rosenbrock4 and adams bashforth moulton1 work", { + skip("Too slow to run") + utils::data(Galapagos_datalist_2types) + pars1 <- c( + 0.195442017, + 0.087959583, + Inf, + 0.002247364, + 0.873605049, + 3755.202241, + 8.909285094, + 14.99999923, + 0.002247364, + 0.873605049, + 0.163 + ) + pars2 <- c(40, 11, 0, 1) + methode <- 'deSolve_R::lsodes' + loglik_lsodes_R <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) + methode <- 'odeint::rosenbrock4' + loglik_rb <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) + methode <- 'odeint::adams_bashforth_moulton_1' + DAISIE_CS_max_steps(100000000) + DAISIE_abm_factor(0.000001) + loglik_abm <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) + expect_equal(loglik_lsodes, loglik_rb) + expect_equal(loglik_lsodes, loglik_abm, tol = 1E-6) +}) From 4a4351584be3116daf9862fc2edf832962385eb6 Mon Sep 17 00:00:00 2001 From: Neves-P Date: Thu, 9 Mar 2023 16:39:33 +0100 Subject: [PATCH 07/13] Fix typo --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 1ab691c1..6363be40 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -110,7 +110,7 @@ Description: Simulates and computes the (maximum) likelihood of a dynamical model of island biota assembly through speciation, immigration and extinction. See e.g. Valente et al. 2015. Ecology Letters 18: 844-852, . -NeedsCompilation: yesc +NeedsCompilation: yes Encoding: UTF-8 VignetteBuilder: knitr URL: https://github.com/rsetienne/DAISIE From 78daa9479cef5b4b7d1a61eeabf026cb6e6372e0 Mon Sep 17 00:00:00 2001 From: Neves-P Date: Thu, 9 Mar 2023 16:47:51 +0100 Subject: [PATCH 08/13] Skip plot tests that would run on headless systems --- tests/testthat/test-DAISIE_plot_sims.R | 5 +++-- tests/testthat/test-DAISIE_plot_stt.R | 1 + 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test-DAISIE_plot_sims.R b/tests/testthat/test-DAISIE_plot_sims.R index c6da7367..93417224 100644 --- a/tests/testthat/test-DAISIE_plot_sims.R +++ b/tests/testthat/test-DAISIE_plot_sims.R @@ -1,5 +1,5 @@ test_that("Example 1", { - + skip("Plots: Run and manually inspect output") data(islands_1type_1000reps) expect_silent( DAISIE_plot_sims( @@ -9,6 +9,7 @@ test_that("Example 1", { }) test_that("Example 2", { + skip("Plots: Run and manually inspect output") data(islands_2types_1000reps) expect_silent( DAISIE_plot_sims( @@ -18,7 +19,7 @@ test_that("Example 2", { }) test_that("Plot plus one", { - + skip("Plots: Run and manually inspect output") data(islands_1type_1000reps) expect_silent( DAISIE_plot_sims( diff --git a/tests/testthat/test-DAISIE_plot_stt.R b/tests/testthat/test-DAISIE_plot_stt.R index b73b1ff0..45a9faea 100644 --- a/tests/testthat/test-DAISIE_plot_stt.R +++ b/tests/testthat/test-DAISIE_plot_stt.R @@ -1,4 +1,5 @@ test_that("use", { + skip("Plots: Run and manually inspect output") utils::data(islands_1type_1000reps, package = "DAISIE") plot_lists <- DAISIE:::DAISIE_convert_to_classic_plot( simulation_outputs = islands_1type_1000reps From 96ee68da0be06bbcdf99a5bff6510c707f1a9473 Mon Sep 17 00:00:00 2001 From: Neves-P Date: Thu, 9 Mar 2023 17:44:13 +0100 Subject: [PATCH 09/13] Add patch notes --- NEWS.md | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/NEWS.md b/NEWS.md index 6effaa8a..bcc74ac3 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,29 @@ # DAISIE 4.3.3 +* Address problem detected by valgrind: unitialized member variable +bulirsch_stoer<>::m_dt_last. + * Patched version of `boost/numeric/odeint/stepper/bulirsch_stoer.hpp`. This + is done by including the patched header file `src/patched_bulrisch_stoer.hpp` + this header before `boost/numeric/odeint` to shadow + `boost/numeriuc/odeint/stepper/bulrisch_stoer.hpp` + The issue is *not* fixed in BOOST_VERSION == 1.81.0. + Must check for fixes in upcomming boost (BH) releases. +* Fix tab spacing in `src/DAISIE_loglik_rhs_FORTRAN.f95`. +* Re-implement clang16 `-D_HAS_AUTO_PTR_ETC` flag fix via Makevars[.win] to +comply with CRAN requests and to force C++ standard C++14 without using +SystemRequirements line in DESCRIPTION, at CRAN's request. +* Skip plotting tests that cause problems in headless systems. These should +be run manually. +* Skip tests of as-of-now experimental steppers from ODEINT rosenbrock4 and +adams bashfort moulton 1 because they are too slow. + +# DAISIE 4.3.2 + +* Apply CRAN suggested fixes to clang16 issues with deprecated C++ functions included the Boost library, which are used in some of the stepper functions. + * Add config.h in a macro, checking for, and setting, `_HAS_AUTO_PTR_ETC` and `BOOST_NO_AUTO_PTR`. + * Change SystemRequirements in DESCRIPTION from C++17 to C++14. +This same fix was applied in package `'DDD'` version 5.2.1. + # DAISIE 4.3.1 * Fix issue that prevented 'covr' from running correctly. From 493942dded9e6e13e34d9bc092d96b56e25a579a Mon Sep 17 00:00:00 2001 From: Neves-P Date: Thu, 9 Mar 2023 17:44:13 +0100 Subject: [PATCH 10/13] patched_bulrisch_stoer.hpp -> patched_bulrisch_stoer.h --- NEWS.md | 24 +++++++++++++++++++ src/DAISIE_odeint.h | 2 +- ...sch_stoer.hpp => patched_bulrisch_stoer.h} | 0 3 files changed, 25 insertions(+), 1 deletion(-) rename src/{patched_bulrisch_stoer.hpp => patched_bulrisch_stoer.h} (100%) diff --git a/NEWS.md b/NEWS.md index 6effaa8a..640f1f99 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,29 @@ # DAISIE 4.3.3 +* Address problem detected by valgrind: unitialized member variable +bulirsch_stoer<>::m_dt_last. + * Patched version of `boost/numeric/odeint/stepper/bulirsch_stoer.hpp`. This + is done by including the patched header file `src/patched_bulrisch_stoer.h` + this header before `boost/numeric/odeint` to shadow + `boost/numeriuc/odeint/stepper/bulrisch_stoer.hpp` + The issue is *not* fixed in BOOST_VERSION == 1.81.0. + Must check for fixes in upcomming boost (BH) releases. +* Fix tab spacing in `src/DAISIE_loglik_rhs_FORTRAN.f95`. +* Re-implement clang16 `-D_HAS_AUTO_PTR_ETC` flag fix via Makevars[.win] to +comply with CRAN requests and to force C++ standard C++14 without using +SystemRequirements line in DESCRIPTION, at CRAN's request. +* Skip plotting tests that cause problems in headless systems. These should +be run manually. +* Skip tests of as-of-now experimental steppers from ODEINT rosenbrock4 and +adams bashfort moulton 1 because they are too slow. + +# DAISIE 4.3.2 + +* Apply CRAN suggested fixes to clang16 issues with deprecated C++ functions included the Boost library, which are used in some of the stepper functions. + * Add config.h in a macro, checking for, and setting, `_HAS_AUTO_PTR_ETC` and `BOOST_NO_AUTO_PTR`. + * Change SystemRequirements in DESCRIPTION from C++17 to C++14. +This same fix was applied in package `'DDD'` version 5.2.1. + # DAISIE 4.3.1 * Fix issue that prevented 'covr' from running correctly. diff --git a/src/DAISIE_odeint.h b/src/DAISIE_odeint.h index ec63abde..a2d531a5 100644 --- a/src/DAISIE_odeint.h +++ b/src/DAISIE_odeint.h @@ -6,7 +6,7 @@ #include "config.h" #include "ublas_types.h" -#include "patched_bulrisch_stoer.hpp" // shadow buggy boost header +#include "patched_bulrisch_stoer.h" // shadow buggy boost header #include #include #include diff --git a/src/patched_bulrisch_stoer.hpp b/src/patched_bulrisch_stoer.h similarity index 100% rename from src/patched_bulrisch_stoer.hpp rename to src/patched_bulrisch_stoer.h From ada6f784022cec4af1842dbb64fa7ff8ba30d309 Mon Sep 17 00:00:00 2001 From: Neves-P Date: Thu, 9 Mar 2023 18:51:35 +0100 Subject: [PATCH 11/13] Fix typo --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index b635dae5..eaf98bc1 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,7 +5,7 @@ bulirsch_stoer<>::m_dt_last. * Patched version of `boost/numeric/odeint/stepper/bulirsch_stoer.hpp`. This is done by including the patched header file `src/patched_bulrisch_stoer.h` before `boost/numeric/odeint` to shadow - `boost/numeriuc/odeint/stepper/bulrisch_stoer.hpp` + `boost/numeric/odeint/stepper/bulrisch_stoer.hpp` The issue is *not* fixed in BOOST_VERSION == 1.81.0. Must check for fixes in upcomming boost (BH) releases. * Fix tab spacing in `src/DAISIE_loglik_rhs_FORTRAN.f95`. From 7121e6dd8c983ec5881f8cdce09cd8de2402e318 Mon Sep 17 00:00:00 2001 From: Neves-P Date: Fri, 10 Mar 2023 11:56:24 +0100 Subject: [PATCH 12/13] Remove test directory and .vscode files --- .vscode/c_cpp_properties.json | 47 ----------------------------------- .vscode/genenv.R | 19 -------------- .vscode/launch.json | 40 ----------------------------- .vscode/tasks.json | 13 ---------- source_odeint.R | 47 ----------------------------------- 5 files changed, 166 deletions(-) delete mode 100644 .vscode/c_cpp_properties.json delete mode 100644 .vscode/genenv.R delete mode 100644 .vscode/launch.json delete mode 100644 .vscode/tasks.json delete mode 100644 source_odeint.R diff --git a/.vscode/c_cpp_properties.json b/.vscode/c_cpp_properties.json deleted file mode 100644 index 19fc441b..00000000 --- a/.vscode/c_cpp_properties.json +++ /dev/null @@ -1,47 +0,0 @@ -// c_cpp_properties.json -// -// C/C++ extension config file -// Some settings to augment IntelliSense for C/C++ files -{ - "configurations": [ - { - "name": "Linux gdb", - "includePath": [ - // recurse into workspace, picks up ./src/ and ./inst/include/ - "${workspaceFolder}/**", - // R headers - "${env:HOME}/opt/bin/Rroot/lib/R/include", - // includes used by our package - "${env:HOME}/opt/bin/Rlibrary/Rcpp/include" - // boost: "${env:HOME}/opt/bin/Rlibrary/BH/include" - ], - "defines": [], - // absolute path to compiler - "compilerPath": "/usr/bin/gcc", - // "compilerPath": "/usr/bin/clang" - "cStandard": "c11", - // should match [[Rcpp::plugins(cpp17)]] - "cppStandard": "c++17" - }, - { - "name": "Linux lldb", - "includePath": [ - // recurse into workspace, picks up ./src/ and ./inst/include/ - "${workspaceFolder}/**", - // R headers - "${env:HOME}/opt/bin/Rroot/lib/R/include", - // includes used by our package - "${env:HOME}/opt/bin/Rlibrary/Rcpp/include" - // boost: "${env:HOME}/opt/bin/Rlibrary/BH/include" - ], - "defines": [], - // absolute path to compiler - "compilerPath": "/usr/bin/clang", - // "compilerPath": "/usr/bin/clang" - "cStandard": "c11", - // should match [[Rcpp::plugins(cpp17)]] - "cppStandard": "c++17" - } - ], - "version": 4 -} diff --git a/.vscode/genenv.R b/.vscode/genenv.R deleted file mode 100644 index 2cd0a945..00000000 --- a/.vscode/genenv.R +++ /dev/null @@ -1,19 +0,0 @@ -# genenv.R -# -# Colects environment variables from R, required by the R *binary*. -# They are normally set by the R *script* which is circumvented if -# we need to run the binary under the debugger as we want to debug R -# not bash... -# -# This script is configured as a task in 'tasks.json' and referenced within -# 'launch.json' as the 'preLaunchTask'. - -env <- Sys.getenv() -envnames <- names(env) -rnames <- envnames[startsWith(envnames, "R_")] -cached_names <- rnames -ld_lib_path <- Sys.getenv("LD_LIBRARY_PATH") -if (ld_lib_path != "") { - cached_names <- c("LD_LIBRARY_PATH", rnames) -} -writeLines(paste0(cached_names, "=", env[cached_names]), ".vscode/.env") diff --git a/.vscode/launch.json b/.vscode/launch.json deleted file mode 100644 index 665df8ae..00000000 --- a/.vscode/launch.json +++ /dev/null @@ -1,40 +0,0 @@ -// launch.json -// -// launch the R *binary* under gdb and run devtools::test() -{ - // Use IntelliSense to learn about possible attributes. - // Hover to view descriptions of existing attributes. - // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387 - "version": "0.2.0", - "configurations": [ - { - "name": "(gbd) devtools::test()", - "type": "cppdbg", - "request": "launch", - // The binary, not the script - "program": "${env:HOME}/opt/bin/Rroot/lib/R/bin/exec/R", - "args": [ - "--vanilla", - "-e", - "devtools::test()" - ], - "stopAtEntry": false, - // needs to be generated, see below - "envFile": "${workspaceFolder}/.vscode/.env", - "cwd": "${workspaceFolder}", - "externalConsole": false, - "MIMode": "gdb", - "setupCommands": [ - { - "description": "Enable pretty-printing for gdb", - "text": "-enable-pretty-printing", - "ignoreFailures": true - } - ], - // 'R' is a script that sets a ton of environment variables - // required by the R binary. This task emulates that part of - // the R script: - "preLaunchTask": "genenv" - } - ] -} diff --git a/.vscode/tasks.json b/.vscode/tasks.json deleted file mode 100644 index e11a2f20..00000000 --- a/.vscode/tasks.json +++ /dev/null @@ -1,13 +0,0 @@ -{ - // See https://go.microsoft.com/fwlink/?LinkId=733558 - // for the documentation about the tasks.json format - "version": "2.0.0", - "tasks": [ - { - "label": "genenv", - "type": "shell", - "command": "Rscript ${workspaceFolder}/.vscode/genenv.R", - "problemMatcher": [] - } - ] -} diff --git a/source_odeint.R b/source_odeint.R deleted file mode 100644 index 5ec7ac01..00000000 --- a/source_odeint.R +++ /dev/null @@ -1,47 +0,0 @@ -# source_odeint.R -# -# partial odeint test suite from tests/testthat/test_odeint.R -# Run with: -# -# R -d "valgrind --tool=memcheck --leak-check=full --error-exitcode=1" --vanilla -e "source('source_odeint.R')" &> source_odeint.log - - -library(DAISIE) - -utils::data(Galapagos_datalist_2types) -pars1 <- c( -0.195442017, -0.087959583, -Inf, -0.002247364, -0.873605049, -3755.202241, -8.909285094, -14.99999923, -0.002247364, -0.873605049, -0.163 -) -pars2 <- c(40, 11, 0, 1) -methode <- 'deSolve_R::lsodes' -loglik_lsodes_R <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) -methode <- 'deSolve_R::lsoda' -loglik_lsoda_R <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) -methode <- 'lsodes' -loglik_lsodes <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) -methode <- 'lsoda' -loglik_lsoda <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) -methode <- 'odeint::runge_kutta_cash_karp54' -loglik_rkck54 <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) -methode <- 'odeint::runge_kutta_fehlberg78' -loglik_rkf78 <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) -methode <- 'odeint::runge_kutta_dopri5' -loglik_rkd5 <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) -methode <- 'odeint::bulirsch_stoer' -loglik_bs <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) -methode <- 'odeint::rosenbrock4' -loglik_rb <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) -methode <- 'odeint::adams_bashforth_moulton_1' -DAISIE_CS_max_steps(100000000) -DAISIE_abm_factor(0.000001) -loglik_abm <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) From b37332753d3b95ab65054b03cb05c4fc7bbb1b59 Mon Sep 17 00:00:00 2001 From: Neves-P Date: Fri, 10 Mar 2023 13:07:59 +0100 Subject: [PATCH 13/13] Update date --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 6363be40..813d99dc 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,7 +2,7 @@ Package: DAISIE Type: Package Title: Dynamical Assembly of Islands by Speciation, Immigration and Extinction Version: 4.3.3 -Date: 2023-02-24 +Date: 2023-03-10 Depends: R (>= 4.1.0) biocViews: Imports: