library(ggplot2)  # package used to visualize results 
library(knitr)    # makes nice tables
library(magrittr) # import the %>% pipeline

library(MEMC)     # MEMC should already be installed, see installation instructions. 

theme_set(theme_bw())  # set a theme to use in all the plots

params <- MEMC::default_params
state  <- MEMC::default_initial
time   <- seq(from = 1, to = 3000, by = 10)

# Set up the model with the LM decay, note this is the default set up and dd.beta is set to a value of 1.
MBdecay_LM <- configure_model(params, state, carbon_pools, carbon_fluxes, MBdecay = "LM", 
                              name = "LM decay")
#> |Model    |DOMdecomp |POMdecomp |MBdecay |
#> |:--------|:---------|:---------|:-------|
#> |LM decay |MM        |MM        |LM      |

# Change the dd.beta parameter (the parameter that controls the strength of the
# density-dependent relationship to a value not equal to 1) 
new_params <- 2
names(new_params) <- "dd.beta"

dd_table <- update_params(new_params, params)
MBdecay_DD <- configure_model(dd_table, state, carbon_pools, carbon_fluxes, MBdecay = "DD", name = "DD decay")
#> |Model    |DOMdecomp |POMdecomp |MBdecay |
#> |:--------|:---------|:---------|:-------|
#> |DD decay |MM        |MM        |DD      |

# Solve the two models
out1 <- solve_model(MBdecay_LM, time)

out2 <- solve_model(MBdecay_DD, time)
out2$param <- 2

out <- rbind(out1, out2, fill = TRUE)

out %>% 
  ggplot(aes(time, value, color = name)) + 
  geom_line() + 
  facet_wrap("variable", scales = "free") + 
  labs(title = "Linear vs Density Dependent Decay", 
       y = unique(out$units), 
       x = "Title") + 
  theme(legend.title = element_blank())

new_params <- 1.5
names(new_params) <- "dd.beta"

new_table  <- update_params(new_params, params)
out3       <- solve_model(MBdecay_DD, time, params = new_table)
out3$param <- 1.5 


new_params <- 2.5
names(new_params) <- "dd.beta"

new_table  <- update_params(new_params, params)
out4       <- solve_model(MBdecay_DD, time, params = new_table)
out4$param <- 2.5

out <- rbind(out2, out3, out4, fill = TRUE)
out$label <- paste(out$name, out$param)

out %>% 
  ggplot(aes(time, value, color = label)) + 
  geom_line() + 
  facet_wrap("variable", scales = "free") + 
  labs(title = "Strength of Density-Dependent MB Decay", 
       y = unique(out$units), 
       x = "Time") + 
  theme(legend.title  = element_blank())