## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = "center",
  fig.width = 7,
  fig.height = 5
)

library(rsofun)
library(dplyr)
library(ggplot2)

## -----------------------------------------------------------------------------
library(rsofun)

# this is to deal with an error p_model_drivers.rds not being found 
p_model_drivers

p_model_validation

## -----------------------------------------------------------------------------
p_model_drivers_vcmax25

p_model_validation_vcmax25

## -----------------------------------------------------------------------------
# Define model parameter values.
# Correspond to maximum a posteriori estimates from Bayesian calibration in
# analysis/02-bayesian-calibration.R.
params_modl <- list(
  kphio              = 5.000000e-02, # chosen to be too high for demonstration
  kphio_par_a        =-2.289344e-03,
  kphio_par_b        = 1.525094e+01,
  soilm_thetastar    = 1.577507e+02,
  soilm_betao        = 1.169702e-04,
  beta_unitcostratio = 146.0,
  rd_to_vcmax        = 0.014,
  tau_acclim         = 20.0,
  kc_jmax            = 0.41
  )

# Run the model for these parameters.
output <- rsofun::runread_pmodel_f(
  p_model_drivers,
  par = params_modl
  )

## -----------------------------------------------------------------------------
# Load libraries for plotting
library(dplyr)
library(tidyr)
library(ggplot2)

# Create data.frame for plotting
df_gpp_plot <- output |>
  tidyr::unnest(data) |>
  dplyr::select(date, gpp_mod = gpp) |>
  dplyr::left_join(
    p_model_validation |>
      tidyr::unnest(data) |>
      dplyr::select(date, gpp_obs = gpp),
    by = "date"
    ) |> 

  # Plot only first year
  dplyr::slice(1:365)

# Plot GPP
ggplot(data = df_gpp_plot) +
  geom_point(
    aes(
      x = date,
      y = gpp_obs,
      color = "Observations"
    ),
  ) +
  geom_line(
    aes(
      x = date,
      y = gpp_mod,
      color = "P-model"
    )
  ) +
  theme_classic() +
  theme(
    panel.grid.major.y = element_line(),
    legend.position = "bottom"
    ) +
  labs(
    x = 'Date',
    y = expression(paste("GPP (g C m"^-2, "s"^-1, ")"))
  ) +
  scale_color_manual(
    NULL,
    breaks = c("Observations",
               "P-model"),
    values = c("black", "tomato"))

## -----------------------------------------------------------------------------
settings <- list(
  method = "GenSA",
  metric = cost_rmse_pmodel,
  control = list(maxit = 100),
  par = list(
    kphio = list(lower=0.02, upper=0.2, init = 0.05)
    )
)

## ----eval=FALSE---------------------------------------------------------------
# # calibrate the model and optimize free parameters
# pars <- calib_sofun(
#   drivers = p_model_drivers,
#   obs = p_model_validation,
#   settings = settings,
#   # extra arguments passed to the cost function:
#   targets = "gpp",             # define target variable GPP
#   par_fixed = params_modl[-1]  # fix non-calibrated parameters to previous
#                                # values, removing kphio
#   )

## ----simulate_calibration_run, include = FALSE--------------------------------
# fake variable as optimization isn't run
pars <- list()
pars$par["kphio"] <- 0.04095957

## -----------------------------------------------------------------------------
# Update the parameter list with calibrated value
params_modl$kphio <- pars$par["kphio"]

# Run the model for these parameters
output_new <- rsofun::runread_pmodel_f(
  p_model_drivers,
  par = params_modl
  )

# Update data.frame for plotting
df_gpp_plot <- df_gpp_plot |> 
  left_join(
    output_new |>
      unnest(data) |>
      select(date, gpp_calib = gpp),
    by = "date"
  )

# Plot GPP
ggplot(data = df_gpp_plot) +
  geom_point(
    aes(
      x = date,
      y = gpp_obs,
      color = "Observations"
    ),
  ) +
  geom_line(
    aes(
      x = date,
      y = gpp_mod,
      color = "P-model (uncalibrated)"
    )
  ) +
  geom_line(
    aes(
      x = date,
      y = gpp_calib,
      color = "P-model (calibrated)"
    )
  ) +
  theme_classic() +
  theme(
    panel.grid.major.y = element_line(),
    legend.position = "bottom"
    ) +
  labs(
    x = 'Date',
    y = expression(paste("GPP (g C m"^-2, "s"^-1, ")"))
  ) +
  scale_color_manual(
    NULL,
    breaks = c("Observations",
               "P-model (uncalibrated)",
               "P-model (calibrated)"),
    values = c("black", "grey", "tomato"))

## -----------------------------------------------------------------------------
run_pmodel_onestep_f_bysite(
  lc4 = FALSE,
  forcing = data.frame(
    temp  = 20,           # temperature, deg C
    vpd   = 1000,         # Pa,
    ppfd  = 300/10^6,     # mol/m2/s
    co2   = 400,          # ppm,
    patm  = 101325        # Pa
  ),
  params_modl = list(
    kphio              = 0.04998,    # setup ORG in Stocker et al. 2020 GMD
    kphio_par_a        = 0.0,        # disable temperature-dependence of kphio
    kphio_par_b        = 1.0,
    beta_unitcostratio = 146.0,
    rd_to_vcmax        = 0.014,      # from Atkin et al. 2015 for C3 herbaceous
    kc_jmax            = 0.41
  ),
  makecheck = TRUE
)
#
#
# make sure to check ?run_pmodel_onestep_f_bysite for units of output quantities.

## -----------------------------------------------------------------------------
library(rpmodel)
out_rpmodel <- rpmodel(
  tc             = 20,           # temperature, deg C
  vpd            = 1000,         # Pa,
  co2            = 400,          # ppm,
  fapar          = 1,            # fraction,
  ppfd           = 300/10^6,     # ~~mol/m2/d~~ or mol/m2/s
                                 # rpmodel docs state that units of ppfd define 
                                 # output units of: lue, gpp, vcmax, rd 
                                 # (as well as vcmax25, gs)
  patm           = 101325,       # Pa
  kphio          = 0.04998,      # quantum yield efficiency as calibrated 
                                 # for setup ORG by Stocker et al. 2020 GMD,
  beta           = 146.0,        # unit cost ratio a/b,
  c4             = FALSE,
  method_jmaxlim = "wang17",
  do_ftemp_kphio = FALSE,        # corresponding to setup ORG
  do_soilmstress = FALSE,        # corresponding to setup ORG
  verbose        = TRUE
  )
  
# out_rpmodel
tidyr::as_tibble(out_rpmodel) |> 
  dplyr::select(vcmax, jmax, vcmax25, jmax25, gs, chi, iwue, rd) |>
  # bring units to same as output of run_pmodel_onestep_f_bysite()
  dplyr::mutate(gs   = gs/1/(300/10^6), # divide by fapar (-) and by ppfd (mol/m2/s)
                iwue = iwue/101325,     # divide by patm (Pa)
                rd   = rd*12            # mutliply by c_molmass (gC/mol)
  )