Merge pull request #158 from rsetienne/develop

v4.4.0

.Rbuildignore

@@ -14,8 +14,8 @@
 
 ^\.Rproj\.user$
 ^\.covrignore$
-^LICENSE$
 ^\.zenodo\.json$
 
 ^\.vscode$
-^source_odeint\.R
+^dummy\.R
+^LICENSE\.md$

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: DAISIE
 Type: Package
 Title: Dynamical Assembly of Islands by Speciation, Immigration and Extinction
-Version: 4.3.4
-Date: 2023-03-16
+Version: 4.4.0
+Date: 2023-03-29
 Depends: R (>= 4.2.0)
 biocViews:
 Imports: 
@@ -105,11 +105,11 @@ Authors@R: c(
            email = "[email protected]",
            role = c("ctb"),
            comment = c(ORCID = "0000-0002-1447-7630")))
-License: GPL-3
+License: GPL (>= 3) | file LICENSE
+Copyright: See the file COPYRIGHTS for various DAISIE copyright details
 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,
-    <DOI:10.1111/ele.12461>.
+    extinction. See Valente et al. (2015) <doi:10.1111/ele.12461>.
 NeedsCompilation: yes
 SystemRequirements: C++17
 Encoding: UTF-8

LICENSE.note

@@ -0,0 +1,9 @@
+The DAISIE package as a whole is distributed under >= GPL-3, the license of which can be found in the distributed file LICENSE. The DAISIE package includes code written by two of the package authors that is distributed under BSL-1.0:
+
+* src/DAISIE_CS.cpp
+* src/DAISIE_IW.cpp
+* src/DAISIE_odeint.h
+* src/DAISIE_types.h
+* src/DAISIE_loglik_rhs_FORTRAN.f95
+
+Full copies of the BSL-1.0 license used by these files is included in `inst/LICENSE_1_0.txt`, as is a license and copyright notice on said files, while details are also in `inst/COPYRIGHTS`.

NAMESPACE

@@ -34,14 +34,14 @@ export(DAISIE_sim_relaxed_rate)
 export(DAISIE_sim_time_dep)
 export(DAISIE_sim_trait_dep)
 export(DAISIE_sim_trait_dep_2K)
-export(are_area_pars)
 export(create_CS_version)
 export(create_area_pars)
 export(create_hyper_pars)
 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%")

NEWS.md

@@ -1,3 +1,27 @@
+# DAISIE 4.4.0
+
+* No longer include patched version of `boost/numeric/odeint/stepper/bulirsch_stoer.hpp`,
+address the issue by passing previously uninitialised variable as `boost::units::quantity<boost::units::si::dimensionless, double>`.
+* Document return type of exported Rcpp functions.
+* Move all internal non-error printing to console `message()` and `warning()`. 
+Add internal functions to address this. `verbose` variable is now numeric, 
+varying from 0 to 3. Increasing values increase amount of messages to be 
+printed. Change default of printing some output to golden rule of silence 
+`verbose == 0`. To print again, set `verbose >= 1`.
+* Remove calls to `options(warn = -1)` that suppress warnings in non standard 
+way. Use `suppressWarnings()` where appropriate.
+* Safely restore graphics settings after plot.
+* Fix references in `DESCRIPTION`.
+* Remove internal (unnecessary) calls to internal functions via `:::`.
+* No longer use roxygen2 tag `internal` to document without index but use
+`noRd` instead. `are_area_pars()` is now internal.
+* Replace `\dontrun` examples in documentation with `\donttest`. Speed up 
+examples.
+* License package C++ source files as BSL-1.0 (c) Hanno Hildenbrandt, FORTRAN 
+files as BSL-1.0 (c) Rampal S. Etienne.
+* Add `LICENSE.note`, `inst/COPYRIGHTS` to clarify license and copyrights. Pipe
+such files in `DESCRIPTION`.
+
 # DAISIE 4.3.4
 
 * Require C++17 via `CXX_STD` flag on Makevars[.win].

R/DAISIE_ML1.R

@@ -61,8 +61,7 @@ DAISIE_loglik_all_choosepar <- function(trparsopt,
       )
     }
     if (is.nan(loglik) || is.na(loglik)) {
-      cat("There are parameter values used
-             which cause numerical problems.\n")
+      message("There are parameter values used which cause numerical problems.")
       loglik <- -Inf
     }
   }
@@ -221,66 +220,56 @@ DAISIE_ML1 <- function(
                                        position = 'lambda_a',
                                        column_to_insert = nc)
   }
-  if (length(namepars[idparsopt]) == 0) {
-    optstr <- "nothing"
-  } else {
-    optstr <- namepars[idparsopt]
-  }
 
-  cat("You are optimizing", optstr, "\n")
-  if (length(namepars[idparsfix]) == 0) {
-    fixstr <- "nothing"
-  } else {
-    fixstr <- namepars[idparsfix]
-  }
-  cat("You are fixing", fixstr, "\n")
+  print_ml_par_settings(
+    namepars = namepars,
+    idparsopt = idparsopt,
+    idparsfix = idparsfix,
+    idparsnoshift = idparsnoshift,
+    all_no_shift = all_no_shift,
+    verbose = verbose
+  )
 
-  if (sum(idparsnoshift %in% (all_no_shift)) != 5) {
-    noshiftstring <- namepars[idparsnoshift]
-    cat("You are not shifting", noshiftstring, "\n")
-  }
   idpars <- sort(c(idparsopt, idparsfix, idparsnoshift, idparseq))
 
   missnumspec <- unlist(lapply(datalist, function(list) {list$missing_species})) # nolint
   if (max(missnumspec) > (res - 1)) {
-    cat(
+    warning(
       "The number of missing species is too large relative to the
-        resolution of the ODE.\n")
+       resolution of the ODE.")
     return(out2err)
   }
 
   if (max(missnumspec) > res/10) {
     warning(
       "The number of missing species is quite low relative to the
-        resolution of the ODE.\n")
+        resolution of the ODE.")
   }
 
   if ((length(idpars) != max(idpars))) {
-    cat("The parameters to be optimized and/or fixed are incoherent.\n")
+    warning("The parameters to be optimized and/or fixed are incoherent.")
     return(out2err)
   }
 
   if ((!all(idpars == 1:max(idpars))) || # nolint
       (length(initparsopt) != length(idparsopt)) ||
       (length(parsfix) != length(idparsfix))) {
-    cat("The parameters to be optimized and/or fixed are incoherent.\n")
+    warning("The parameters to be optimized and/or fixed are incoherent.")
     return(out2err)
   }
   if (length(idparseq) == 0) {
   } else {
     if (ddmodel == 3) {
-      cat("Equilibrium optimization is not implemented for ddmodel = 3\n")
+      warning("Equilibrium optimization is not implemented for ddmodel = 3")
     } else {
-      cat(
+      message(
         "You are assuming equilibrium. Extinction and/or immigration will
           be considered a function of the other parameters, the species
           pool size, the number of endemics,
-          and/or the number of non-endemics\n"
+          and/or the number of non-endemics"
       )
     }
   }
-  cat("Calculating the likelihood for the initial parameters.", "\n")
-  utils::flush.console()
   trparsopt <- initparsopt / (1 + initparsopt)
   trparsopt[which(initparsopt == Inf)] <- 1
   trparsfix <- parsfix / (1 + parsfix)
@@ -313,20 +302,17 @@ DAISIE_ML1 <- function(
     abstolint = tolint[1],
     reltolint = tolint[2]
   )
-  cat(
-    "The loglikelihood for the initial parameter values is",
-    initloglik,
-    "\n"
-  )
+
+  print_init_ll(initloglik = initloglik, verbose = verbose)
+
   if (initloglik == -Inf) {
-    cat(
+    warning(
       "The initial parameter values have a likelihood that is equal to 0 or
-       below machine precision. Try again with different initial values.\n"
+       below machine precision. Try again with different initial values."
     )
     return(out2err)
   }
-  cat("Optimizing the likelihood - this may take a while.", "\n")
-  utils::flush.console()
+
   out <- DDD::optimizer(
     optimmethod = optimmethod,
     optimpars = optimpars,
@@ -347,9 +333,9 @@ DAISIE_ML1 <- function(
     num_cycles = num_cycles
   )
   if (out$conv != 0) {
-    cat(
+    warning(
       "Optimization has not converged.
-        Try again with different initial values.\n")
+       Try again with different initial values.")
     out2 <- out2err
     out2$conv <- out$conv
     return(out2)
@@ -414,7 +400,7 @@ DAISIE_ML1 <- function(
       conv = unlist(out$conv)
     )
     pars_to_print <- MLpars1[1:6]
-    parnames <- c('lambda^c','mu','K','gamma','lambda^a','prob_init_pres')
+    parnames <- c("lambda^c", "mu", "K", "gamma", "lambda^a", "prob_init_pres")
   } else {
     out2 <- data.frame(
       lambda_c = MLpars1[1],
@@ -431,17 +417,16 @@ DAISIE_ML1 <- function(
   }
   print_parameters_and_loglik(pars = pars_to_print,
                               loglik = ML,
-                              verbose = TRUE,
+                              verbose = verbose,
                               parnames = parnames,
                               type = 'island_ML')
   if (eqmodel > 0) {
     M <- calcMN(datalist, MLpars1)
     ExpEIN <- DAISIE_ExpEIN(datalist[[1]]$island_age, MLpars1, M) # nolint start
-    cat("The expected number of endemics, non-endemics, and the total at
-        these parameters is: ",
-        ExpEIN[[1]],
-        ExpEIN[[2]],
-        ExpEIN[[3]]
+    message(
+      paste0("The expected number of endemics, non-endemics, and the total at ",
+        "these parameters is: "),
+      paste(ExpEIN[[1]], ExpEIN[[2]], ExpEIN[[3]])
     ) # nolint end
   }
   return(invisible(out2))

R/DAISIE_ML2.R

@@ -30,7 +30,7 @@ DAISIE_loglik_all_choosepar2 <- function(
       }
    }
    if (is.nan(loglik) || is.na(loglik)) {
-      cat("There are parameter values used which cause numerical problems.\n")
+      warning("There are parameter values used which cause numerical problems.")
       loglik <- -Inf
    }
    return(loglik)
@@ -120,8 +120,6 @@ DAISIE_ML2 <- function(
 # - cond = conditioning; currently only cond = 0 is possible
 #  . cond == 0 : no conditioning
 #  . cond == 1 : conditioning on presence on the island
-
-  options(warn = -1)
   out2err <- data.frame(lambda_c = NA, mu = NA, K = NA, gamma = NA, lambda_a = NA, loglik = NA, df = NA, conv = NA)
   out2err <- invisible(out2err)
   numisl <- length(datalist)
@@ -131,31 +129,30 @@ DAISIE_ML2 <- function(
   }
 
   if (missnumspec > (res - 1)) {
-    cat("The number of missing species is too large relative to the resolution of the ODE.\n")
+    warning("The number of missing species is too large relative to the resolution of the ODE.")
     return(out2err)
   }
   if (all((sort(unique(as.vector(idparsmat))) != sort(c(idparsopt, idparsfix)))) ||
       (length(initparsopt) != length(idparsopt)) ||
       (length(parsfix) != length(idparsfix))) {
-    cat("The parameters to be optimized and/or fixed are incoherent.\n")
+    warning("The parameters to be optimized and/or fixed are incoherent.")
     return(out2err)
   }
-  cat("Calculating the likelihood for the initial parameters.", "\n")
-  utils::flush.console()
+
   trparsopt <- initparsopt / (1 + initparsopt)
   trparsopt[which(initparsopt == Inf)] <- 1
   trparsfix <- parsfix / (1 + parsfix)
   trparsfix[which(parsfix == Inf)] <- 1
   pars2 <- c(res, ddmodel, cond, 0, island_ontogeny)
   optimpars <- c(tol, maxiter)
   initloglik <- DAISIE_loglik_all_choosepar2(trparsopt = trparsopt, trparsfix = trparsfix, idparsopt = idparsopt, idparsfix = idparsfix, idparsmat = idparsmat, pars2 = pars2, datalist = datalist, methode, abstolint = tolint[1], reltolint = tolint[2])
-  cat("The loglikelihood for the initial parameter values is", initloglik, "\n")
+
+  print_init_ll(initloglik = initloglik, verbose = verbose)
+
   if (initloglik == -Inf) {
-    cat("The initial parameter values have a likelihood that is equal to 0 or below machine precision. Try again with different initial values.\n")
+    warning("The initial parameter values have a likelihood that is equal to 0 or below machine precision. Try again with different initial values.")
     return(out2err)
   }
-  cat("Optimizing the likelihood - this may take a while.", "\n")
-  utils::flush.console()
   out <- DDD::optimizer(
     optimmethod = optimmethod,
     optimpars = optimpars,
@@ -174,7 +171,7 @@ DAISIE_ML2 <- function(
     num_cycles = num_cycles
   )
   if (out$conv != 0) {
-    cat("Optimization has not converged. Try again with different initial values.\n")
+    warning("Optimization has not converged. Try again with different initial values.")
     out2 <- out2err
     out2$conv <- out$conv
     return(out2err)
@@ -206,7 +203,7 @@ DAISIE_ML2 <- function(
                      conv = unlist(out$conv))
   print_parameters_and_loglik(pars = MLpars1,
                               loglik = ML,
-                              verbose = TRUE,
+                              verbose = verbose,
                               type = 'multiple_island_ML')
   return(invisible(out2))
 }