merge symmetry project code (#39)

* Update bibliography and finalizing manuscript.

* Final draft

* Add event_definition = NULL

* Integrate project from @pascalbartschi

Functions were extended with the 'symmetric' parameter as well as a new ODI function and two mote plot functions.

* Delete manuscript_bibliography.bib

* v0.9.1.5 update

* v0.9.1.1: Changes from Pascal

* v0.9.1.1: Symmetric

* new_runsim_parameter list now holds sym_axis value

* added function

* added run_simulation_symmetric check and documentation

* added forcing of aO AND aS, as well as updated documentation

* event definiton needed for syymmetric simulation. 2 log10a forcing functions

* adapted for change of aO and aS

* newline at end of file

* newline

* - runs ssfind simulation compatible with symmetric bushplus simulation
- updated documentation: this function is used ti find ss, but doesn't find per se!

* function to plot result of ssfind simulation (sym and other)

* renamed function

* inserted compbility for symmetry simulation

* updated docucomentation

* get symmmetry measures function

* added documentation of get sym measures

* plotting function

* functionand documnentation added

* upated importFrom

* convience function to plot symmetric trajectories

* update for stacking returned plots: use of ggpubr::arrange as '/' was depricated in R 4.2.3

* removed geom_vline of sym_axis

* developed plotting functions for higher accuracy and broader universality in plotting magnitudes

* developed plotting functions

* developed package for simulation with asymmetric stressor pressure on CB/SB, controlled by parameter asym_factor

* bug fixes: asymmetric stressor manipulation

* secondary axis plotting implemented

* broader functionality, bug fixes, preset args

* bug fixes

* plot label changed

* removed ggarrange, reintroduced patchwork

* bug fix asymmetryic 'log10a_series' initialization

* bug fix secondary x-axis

* x label fix

* switched to colorblind compatible hex colors

* refined mirroring of diffusivity vectors

* started symmmetric User-Guide

* renamed function to set_diffusivities

* bugfix in if-condition of set_diffusivties.R

* added data to inst/data-pkg and devoloped sym-user-guide to have slider

* renamed sym measurememnts functions

* finished first draft of symsys-user-guide

* added .RDS file for sym-sys user guide

* bugfix in variable name inconsistency

* bugfix function call plot_trajectory_symmetry.R

* plot labels have units

* bug fix

* vignettes working and checked

* update version number to 0.9.2

* bug fix v0.9.2

* small changes

* updated sym user guide, bugfix in set_diffusivities

* removed hardcoding of theme_bw()

* final formatting of SymSys-User-Guide

* Finalises for merge
Fixes:
 - missing global variable definitiuons
 - missing global functions in import
 - consistent formating
 - other non-functioning related edits

---------

Co-authored-by: Rainer M Krug <Rainer@krugs.de>
Co-authored-by: Owen Petchey <owen.petchey@ieu.uzh.ch>

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: microxanox
-Version: 0.9.1
-Date: 2022-09-05
-Title: Oxic-anoxic regime shifts in microbial communities
+Version: 0.9.3
+Date: 2023-11-13
+Title: Oxic-Anoxic Regime Shifts in Microbial Communities
 Description: Model to simulate a three functional group system with four chemical substrates
   using a set of ordinary differential equations. 
   Simulations can be run individually or over a parameter range, to find stable states. The 
@@ -12,15 +12,20 @@ Authors@R:
     c(person(given = "Owen L.",
              family = "Petchey",
              role = c("aut", "cre"),
-             email = "Owen.Petchey@uzh.ch",
+             email = "Owen.Petchey@ieu.uzh.ch",
              comment = c(ORCID = "0000-0002-7724-1633")),
       person(given = "Rainer M.",
              family = "Krug",
              role = "ctb",
              email = "Rainer.Krug@uzh.ch",
-             comment = c(ORCID = "0000-0002-7490-0066"))
+             comment = c(ORCID = "0000-0002-7490-0066")),
+      person(given = "Pascal",
+             family = "Baertschi",
+             role = "ctb",
+             email = "compbio.baertschi@gmail.com",
+             comment = c(ORCID = "0000-0002-6533-0399"))
       )
-Maintainer: Owen L Petchey <owen.petchey@ieu.uzh.ch>
+Maintainer: Owen L. Petchey <Owen.Petchey@ieu.uzh.ch>
 URL: https://github.com/UZH-PEG/microxanox/
 BugReports: https://github.com/UZH-PEG/microxanox/issues 
 License: MIT + file LICENSE
@@ -38,14 +43,21 @@ Imports:
   dplyr,
   tidyr,
   stringr,
-  multidplyr
+  multidplyr,
+  ggpubr,
+  DescTools
 Suggests: 
     knitr,
     rmarkdown,
     rootSolve,
     tidyverse,
-    testthat (>= 3.0.0)
-RoxygenNote: 7.2.1
+    testthat (>= 3.0.0),
+    roxyglobals
+RoxygenNote: 7.2.3
 Encoding: UTF-8
 VignetteBuilder: knitr
 Config/testthat/edition: 3
+Roxygen: list(roclets = c("collate", "namespace", "rd",
+    "roxyglobals::global_roclet"))
+Config/roxyglobals/filename: globals.R
+Config/roxyglobals/unique: FALSE

Makefile

@@ -90,13 +90,13 @@ deps:
 
 ####
 
-_docs:
-	Rscript -e "devtools::document(roclets = c('rd', 'collate', 'namespace', 'vignette', 'roxyglobals::global_roclet'))"
+doc:
+	Rscript -e "devtools::document()"
 
-_codemeta:
+metadata:
 	Rscript -e "codemetar::write_codemeta()"
 
-docs: _docs _codemeta 
+docs: docs metadata 
 
 ####
 
@@ -111,46 +111,34 @@ clean_vignettes:
 ####
 
 
-_build: 
-	cd ..;\
-	R CMD build $(PKGSRC)
-
-build: docs _build
+build: doc
 	cd ..;\
 	R CMD build $(PKGSRC)
 
 ####
 
-_drat:
+drat:
 	cd
 	@Rscript -e "drat::insertPackage('./../$(PKGNAME)_$(PKGVERS).tar.gz', repodir = './../../drat/', commit = TRUE)"
-
-drat: docs _build _drat
-
+
 ####
 
-_build-cran:
+build-cran: doc
 	cd ..;\
 	R CMD build $(PKGSRC)
 
-build-cran: docs _build_cran
-
 ####
 
-_install:
+install:
 	cd ..;\
 	R CMD INSTALL $(PKGNAME)_$(PKGVERS).tar.gz
 
-install: build _install
-
 ####
 
-_check:
+check: build-cran
 	cd ..;\
 	R CMD check $(PKGNAME)_$(PKGVERS).tar.gz --as-cran
 
-check: build-cran _check
-
 clean_check:
 	$(RM) -r ./../$(PKGNAME).Rcheck/
 

NAMESPACE

@@ -6,6 +6,7 @@ export(add_strain_var)
 export(bushplus_dynamic_model)
 export(event_definition_1)
 export(event_definition_2)
+export(event_definition_symmetric)
 export(get_final_states_a_N)
 export(get_hysteresis_max)
 export(get_hysteresis_min)
@@ -16,6 +17,7 @@ export(get_ssbyaN_parameter)
 export(get_stability_measures)
 export(get_stability_measures_replication_ssfind_result)
 export(get_stability_measures_temporal_ssfind_result)
+export(get_symmetry_measurements)
 export(growth1)
 export(growth2)
 export(inhibition)
@@ -30,34 +32,59 @@ export(new_runsim_results)
 export(new_strain_parameter)
 export(new_temporal_ssfind_results)
 export(plot_dynamics)
+export(plot_dynamics_symmetric)
+export(plot_symmetry_measures)
+export(plot_temporal_ss)
+export(plot_trajectory_symmetry)
 export(run_replication_ssfind)
 export(run_simulation)
+export(run_simulation_symmetric)
 export(run_temporal_ssfind)
 export(run_temporal_ssfind_experiment)
+export(run_temporal_ssfind_symmetric)
+export(set_diffusivities)
 export(set_temporal_ssfind_initial_state)
+export(symmetric_bushplus_dynamic_model)
 import(patchwork)
+importFrom(DescTools,AUC)
 importFrom(deSolve,ode)
 importFrom(dplyr,across)
+importFrom(dplyr,all_of)
 importFrom(dplyr,case_when)
 importFrom(dplyr,collect)
 importFrom(dplyr,filter)
 importFrom(dplyr,mutate)
 importFrom(dplyr,pull)
+importFrom(dplyr,recode)
 importFrom(dplyr,rename)
 importFrom(dplyr,row_number)
 importFrom(dplyr,rowwise)
 importFrom(dplyr,select)
 importFrom(dplyr,slice)
+importFrom(dplyr,starts_with)
 importFrom(dplyr,summarise)
 importFrom(ggplot2,aes)
+importFrom(ggplot2,arrow)
 importFrom(ggplot2,facet_wrap)
 importFrom(ggplot2,geom_line)
+importFrom(ggplot2,geom_path)
+importFrom(ggplot2,geom_point)
+importFrom(ggplot2,geom_ribbon)
+importFrom(ggplot2,geom_segment)
 importFrom(ggplot2,ggplot)
 importFrom(ggplot2,guide_legend)
 importFrom(ggplot2,guides)
+importFrom(ggplot2,labs)
+importFrom(ggplot2,scale_color_manual)
 importFrom(ggplot2,scale_colour_manual)
+importFrom(ggplot2,scale_fill_manual)
+importFrom(ggplot2,scale_x_continuous)
+importFrom(ggplot2,sec_axis)
+importFrom(ggplot2,unit)
 importFrom(ggplot2,xlab)
 importFrom(ggplot2,ylab)
+importFrom(ggplot2,ylim)
+importFrom(ggpubr,ggarrange)
 importFrom(grDevices,colorRampPalette)
 importFrom(magrittr,"%>%")
 importFrom(multidplyr,cluster_library)
@@ -66,6 +93,7 @@ importFrom(multidplyr,partition)
 importFrom(stats,approx)
 importFrom(stats,approxfun)
 importFrom(stats,lm)
+importFrom(stats,median)
 importFrom(stats,na.omit)
 importFrom(stats,predict)
 importFrom(stringr,str_detect)

NEWS.md

@@ -1,5 +1,7 @@
+# 0.9.2: Added symmetric model functionality
+
 # 0.9.1: Fix typo in bushplus_dynamic_model.R
-The error had no cincernible impact on the results of the vignettes
+The error had no discernible impact on the results of the vignettes
 
 # 0.9.0: For first submission submission
 

R/0.microxanox-package.R

@@ -16,7 +16,7 @@
 #' - A user guide.
 #' - A reproduction of some results of Bush et al 2017, on which the simulation is based. It shows the alternate stable states. It also extends the original study by adding temporal forcing and therefore temporal switching between states, and also effects of biotic composition on the response to environmental change.
 #'
-#' 
+#'
 #' @name microxanox-package
 #' @aliases microxanox
 #' @docType package

R/add_strain_var.R

@@ -26,6 +26,9 @@
 #'   mentioned variables. 
 #'
 #' @return the value of strain_parameter with added strain variability
+#' 
+#' @autoglobal
+#' 
 #' @export
 #'
 #' @md

R/bushplus_dynamic_model.R

@@ -1,5 +1,5 @@
 #' The rate equations, as published in the Bush et al 2017
-#' \url{https://doi.org/10.1093/clinchem/39.5.766} paper, but with forcing of
+#' \doi{10.1093/clinchem/39.5.766} paper, but with forcing of
 #' oxygen diffusivity `a_0` potential added, and the possibility to simulate
 #' multiple strains per functional group
 #'
@@ -10,20 +10,20 @@
 #'   `new_runsim_parameter()``
 #' @param log10a_forcing_func function to change oxygen diffusivity `a` depending on `t`
 #' @param ... not used. Needed to catch additional parameters.
-#' 
+#'
 #' @return a list containing two elements, namely the rate of change of the
 #'   strains, and also the current values of oxygen diffusivity `a`.
 #' @md
-#' 
+#'
+#' @autoglobal
+#'
 #' @export
 bushplus_dynamic_model <- function(
-  t, 
-  state, 
-  parameters, 
-  log10a_forcing_func,
-  ...
-){
-
+    t,
+    state,
+    parameters,
+    log10a_forcing_func,
+    ...) {
   ## unpack state variables from the state object
   CB <- state[grep("CB", names(state))]
   names_CB <- names(CB)[order(names(CB))]
@@ -35,8 +35,8 @@ bushplus_dynamic_model <- function(
   names_SB <- names(SB)[order(names(SB))]
   SB <- as.numeric(SB[order(names(SB))])
   # print(c(CB, PB, SB))
-  
-  
+
+
   # CB rates of change
   CB_growth_rate <- growth1(state["P"], parameters$CB$g_max_CB, parameters$CB$k_CB_P) * inhibition(state["SR"], parameters$CB$h_SR_CB) * CB
   CB_mortality_rate <- parameters$CB$m_CB * CB
@@ -46,32 +46,32 @@ bushplus_dynamic_model <- function(
   PB_growth_rate <- growth2(state["P"], state["SR"], parameters$PB$g_max_PB, parameters$PB$k_PB_P, parameters$PB$k_PB_SR) * inhibition(state["O"], parameters$PB$h_O_PB) * PB
   PB_mortality_rate <- parameters$PB$m_PB * PB
   PB_rate <- PB_growth_rate - PB_mortality_rate + parameters$PB$i_PB
-  
+
   # SB rates of change
   SB_growth_rate <- growth2(state["P"], state["SO"], parameters$SB$g_max_SB, parameters$SB$k_SB_P, parameters$SB$k_SB_SO) * inhibition(state["O"], parameters$SB$h_O_SB) * SB
   SB_mortality_rate <- parameters$SB$m_SB * SB
   SB_rate <- SB_growth_rate - SB_mortality_rate + parameters$SB$i_SB
-  
+
   # Substrate rates of change
   SO_rate <- sum(1 / parameters$PB$y_SR_PB * PB_growth_rate) -
     sum(1 / parameters$SB$y_SO_SB * SB_growth_rate) +
     parameters$c * state["O"] * state["SR"] +
-    parameters$a_S * (parameters$back_SO -state[["SO"]])
-  
-  SR_rate <- - sum(1 / parameters$PB$y_SR_PB * PB_growth_rate) +
+    parameters$a_S * (parameters$back_SO - state[["SO"]])
+
+  SR_rate <- -sum(1 / parameters$PB$y_SR_PB * PB_growth_rate) +
     sum(1 / parameters$SB$y_SO_SB * SB_growth_rate) -
     parameters$c * state["O"] * state["SR"] +
     parameters$a_S * (parameters$back_SR - state["SR"])
-  
+
   O_rate <- sum(parameters$CB$Pr_CB * CB_growth_rate) -
     parameters$c * state["O"] * state["SR"] +
     10^log10a_forcing_func(t) * (parameters$back_O - state["O"])
-  
-  P_rate <- - sum(1 / parameters$CB$y_P_CB * CB_growth_rate) -
+
+  P_rate <- -sum(1 / parameters$CB$y_P_CB * CB_growth_rate) -
     sum(1 / parameters$PB$y_P_PB * PB_growth_rate) -
     sum(1 / parameters$SB$y_P_SB * SB_growth_rate) +
     parameters$a_P * (parameters$back_P - state["P"])
-  
+
   # Assemble results
   result <- list(
     c(
@@ -85,7 +85,7 @@ bushplus_dynamic_model <- function(
     ),
     a = log10a_forcing_func(t)
   )
-  
+
   # Name results
   names(result[[1]]) <- c(
     parameters$CB$strain_name,
@@ -96,7 +96,6 @@ bushplus_dynamic_model <- function(
     "O_rate",
     "P_rate"
   )
-  
+
   return(result)
 }
-